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ABSTRACT 


This  thesis  is  concerned  with  the  control  of  batch  reactions  which 
contain  uncertain  parameters.  Simulation  results  are  presented  for  two 
reaction  systems,  a  single  reversible  reaction  and  the  complex  reaction 
system  of  Rosenbrock  and  Storey  [10]. 

Optimal  and  suboptimal  control  policies  were  studied  and  compared 
under  ideal  conditions  of  noise  free  measurements  and  known  kinetic  para¬ 
meters.  The  suboptimal  control  policies  resulted  in  a  final  yield  close 
to  the  values  for  the  optimal  control  policy.  The  parameter  sensitivity 
of  optimal  and  suboptimal  control  policies  was  also  investigated  to 
determine  the  feasibility  of  the  various  control  policies  when  parameters 
are  changed  in  an  unanticipated  manner. 

To  cope  with  the  problem  of  unknown  parameters  or  undetected  para¬ 
meter  variations,  a  nonlinear  filter  technique  was  used  to  perform  on¬ 
line  parameter  estimation.  To  reduce  the  computational  requirements 
associated  with  on-line  parameter  estimation,  two  simplified  nonlinear 
filter  algorithms  were  also  investigated.  It  was  found  that  the  simpli¬ 
fied  algorithms  provide  a  feasible  approach  for  parameter  estimation 
and  resulted  in  significant  reduction  in  computation  time.  However,  the 
initial  covariance  matrix  £(0)  must  usually  be  determined  by  trial  and 
error. 

The  robustness  of  the  proposed  nonlinear  filter  algorithms  was 
verified  by  investigating  the  effect  of  data  sampling  interval,  noise 
level  and  inaccurate  initial  estimates  of  states.  The  simplified  filter 
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algorithms  were  at  least  as  robust  as  the  original  filter. 

The  use  of  a  nonlinear  filter  in  the  feedback  control  resulted  in 
improved  control  when  the  kinetic  parameters  changed  and  were  "tracked" 
by  the  filter. 
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CHAPTER  I 
INTRODUCTION 

Although  continuous  operation  offers  many  well-known  advantages 
many  industrial  processes  are  still  run  on  a  batch  basis.  Conditions 
which  favor  batch  operation  include  low  volume,  multiple  step  pro¬ 
cessing,  low  rates  of  reaction,  and  the  desire  for  flexibility  when 
product  obsolescence  is  of  concern  [1].  In  the  chemical  industry, 
batch  reactors  are  widely  used  in  the  manufacture  of  fine  chemicals, 
pharmaceuticals  and  polymers. 

The  control  and  automation  problems  associated  with  batch 
reactors  differ  significantly  from  those  of  continuous  reactors  in 
several  ways  [2]: 

(1)  Batch  operation  is  inherently  unsteady  state  with 
reactor  conditions  changing  over  a  wide  range  of 
values.  Reactor  start  up  is  a  daily  routine  rather 
than  an  infrequent  occurrence,  as  in  continuous 
processes . 

(2)  A  batch  process  typically  consists  of  a  series  of 
sequential  steps  with  relatively  large  labor  require 
ments.  Consequently,  automation  will  require  a 
logical  sequencing  system  and  a  digital  computer  is 
an  obvious  candidate  if  the  necessary  economic 
justification  can  be  made. 

(3)  Physical  and  analytical  measurements  are  much  more 
difficult  to  obtain  in  batch  processes.  Because 
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reactor  conditions  vary  widely  during  batch  opera¬ 
tion,  the  measurement  range  is  large  and  yet  accuracy 
requirements  for  measurements  are  higher  in  batch 
processes.  By  contrast,  for  a  well -tuned  continuous 
process  inaccurate  measurements  can  be  tolerated  as 
long  as  a  reasonable  degree  of  reproducibility  is 
attained  [1]. 

For  continuous  processes,  well-known  control  system  design 
techniques  can  be  employed  providing  that  a  satisfactory  linear  (or 
linearized)  process  model  is  available.  In  particular,  single¬ 
variable  feedback  control  systems  based  on  transfer  function  process 
models  have  been  widely  used  in  the  process  industries.  However,  for 
batch  reactors  the  inherently  nonlinear  chemical  kinetics  and  wide 
range  of  reactor  operating  conditions  discourage  the  use  of  a  linear 
process  model  such  as  a  transfer  function. 

An  important  consideration  in  designing  a  batch  reactor  con¬ 
trol  system  is  the  sensitivity  of  the  control  system  performance  to 
the  assumed  values  of  the  process  parameters.  Reaction  rate  expres¬ 
sions  which  are  accurate  over  a  wide  range  of  conditions  are  very 
difficult  to  determine.  Furthermore,  for  batch  reactors,  the  char¬ 
acteristics  of  the  initial  charge  may  change  unpredictably  from  batch 
to  batch  due  to  variations  in  feed  composition,  catalyst  properties, 
or  improper  cleaning  of  the  reactor  after  the  previous  cycle  of 
operation.  The  central  concern  of  this  study  is  to  determine  a 
feasible  control  policy  for  batch  reactors  when  process  parameters 
change  unpredictably  from  batch  to  batch. 
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The  optimal  control  of  chemical  reactors  has  been  studied 
extensively  during  the  last  decade  with  definitive  contributions  being 
made  by  Aris  and  coworkers  [3,4,5 ,6,7 ,8,9] .  In  this  approach  the 
objective  is  to  determine  the  control  policy  which  will  maximize  (or 
minimize)  a  specific  "performance  index".  Thus  a  dynamic  optimization 
problem  is  formulated  and  the  solution  obtained  by  utilization  of 
optimization  techniques  such  as  dynamic  programming,  Pontryagin's 
maximum  principle,  or  hi  1 1 -cl imbing  techniques  [10],  For  batch 
reactors  an  obvious  performance  index  is  the  minimum  time  required  to 
achieve  a  certain  product  yield  or  equivalently  [3,11],  the  maximum 
yield  which  can  be  obtained  in  a  specific  time  interval. 

However,  optimal  control  policies  for  batch  reactors  contain 
several  disadvantages.  Except  for  special  cases  such  as  a  single 
reversible  reaction,  the  optimal  control  law  is  typically  "open  loop" 
in  that  the  manipulated  variables  are  specified  as  a  function  of  time 
rather  than  as  a  function  of  the  measured  process  variables,  as  in  a 
feedback  control  law.  Consequently,  process  disturbances  or 
unanticipated  changes  in  process  parameters  may  have  a  detrimental 
effect  on  control  system  performance.  As  noted  above,  for  batch 
reactions  parameter  values  may  be  uncertain  or  change  unpredi ctably 
from  batch  to  batch.  Even  if  process  parameters  could  be  estimated 
"on-line"  (as  proposed  below),  updating  the  optimal  control  policy 
based  on  current  parameter  values  is  not  attractive  on  account  of 
the  extensive  calculations  required  to  determine  the  optimal  control 
policy  for  a  nonlinear  system.  Such  calculations  would  not  be 
practical  for  most  process  control  computers. 
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Alternatively,  several  "suboptimal"  control  schemes  incor¬ 
porating  feedback  control  action  have  been  proposed  for  batch 
reactors  [12,13,14,15]=  Millman  and  Katz  [13]  have  studied  the  con¬ 
ventional  three  mode  controller  (proportional -integral -deri vati ve) 
for  several  types  of  simulated  consecutive  and  parallel  reaction 
systems.  Their  examples  indicated  that  the  use  of  controller  con¬ 
stants  determined  by  hill-climbing  techniques  resulted  in  reaction 
yields  which  were  very  close  to  the  optimal  yields.  However,  their 
suboptimal  control  policy  would  also  be  difficult  to  adapt  on-line 
to  accomodate  changing  conditions  because  the  required  calculations 
are  quite  extensive. 

A  second  class  of  suboptimal  control  policies  for  nonlinear 
systems  is  concerned  with  minimizing  a  performance  index  based  on 
instantaneous  process  conditions  rather  than  the  process  conditions 
during  the  entire  time  period.  The  resulting  control  scheme  usually 
includes  feedback  action  and  can  easily  be  updated  if  process  para¬ 
meter  changes  are  detected  by  some  method  of  parameter  estimation. 
Suboptimal  control  schemes  based  on  minimizing  an  instantaneous 
performance  index  have  been  investigated  by  Seinfeld  both  for  contin¬ 
uous  reactors  [16]  and  batch  reactors  [11]. 

In  this  study  the  problem  of  interest  is  the  determination  of 
control  policies  for  batch  reactors  whose  process  parameters  change 
in  an  unpredictable  fashion  from  batch  to  batch.  The  general  approach 
adopted  is  to  assume  a  process  model  of  a  particular  form  and  esti¬ 
mate  model  parameters  on-line  from  noisy  measurements  using 
sequential  estimation  techniques  (i.e.  nonlinear  filtering).  The 
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updated  parameter  estimates  are  continuously  incorporated  into  the 
control  law  and  thus  estimation  and  control  proceed  simultaneously. 

This  approach  is  often  referred  to  as  an  adaptive  control  scheme 
[5, 6, 7, 8].  A  significant  advantage  of  this  approach  is  that  esti¬ 
mates  of  inaccessible  state  variables  such  as  chemical  composition 
are  automatically  generated  during  parameter  estimation.  This  per¬ 
mits  the  employment  of  "inferential  control"  techniques  which  use 
estimated  values  of  unmeasured  state  variables  for  control.  Such 
techniques  have  been  successfully  implemented  by  Fisher  and 
Jacobson  [17]  for  a  pilot  plant  evaporator. 

A  major  disadvantage  of  the  proposed  approach  (i.e.  on-line 
parameter  estimation  plus  control)  is  that  implementation  will 
inevitably  require  the  availability  of  a  process  control  computer. 
However,  successful  industrial  applications  of  computer  control 
techniques  for  batch  reactors  have  recently  been  reported  [2,18,19,20]. 
The  main  justification  for  computer  control  of  batch  processes  is  the 
resulting  reduced  labor  costs  and  improved  productivity  [1].  Other 
advantages  include  improved  process  reproducibility  and  precise 
documentation  [21]. 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter 
II  presents  a  discussion  of  existing  nonlinear  filtering  techniques 
with  emphasis  on  simplified  filters  with  reduced  computational 
requirements.  The  combined  parameter  estimation  and  control  scheme 
is  applied  to  a  single  reaction  system  in  Chapter  III  and  to  a  complex 
reaction  system  in  Chapter  IV.  The  parameter  sensitivity  of  both 
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optimal  and  suboptimal  control  schemes  is  investigated  and  the  robust¬ 
ness  of  the  proposed  scheme  to  noise  level,  data  sampling  intervals 
and  other  parameters  is  also  considered.  Finally,  the  conclusions 
of  this  study  are  presented  in  Chapter  V. 


7 


CHAPTER  II 

NONLINEAR  FILTERING  TECHNIQUES 

2.1  Introduction 

Mathematical  models  of  batch  reactors  typically  consist  of 
material  and  energy  balances  written  in  the  form  of  a  set  of  nonlinear 
ordinary  differential  equations.  These  models  are  far  from  precise 
due  to  the  uncertainty  associated  with  the  values  of  process  para¬ 
meters  such  as  rate  constants  and  heat  transfer  coefficients,  and 
variations  in  the  characteristics  of  the  initial  charge  from  batch 
to  batch.  In  this  analysis  it  will  be  assumed  that  the  effects  of 
such  modelling  errors  can  be  significantly  reduced  by  considering 
certain  parameters  in  the  set  of  ordinary  differential  equations  to 
be  "unknown".  The  task  of  the  adaptive  control  scheme  is  then  to 
estimate  the  unknown  parameters  from  noisy  process  measurements 
during  a  batch  run,  and  to  make  the  appropriate  corrections  in  con¬ 
trol  action  to  compensate  for  changes  in  process  parameters. 

The  general  problem  of  estimating  parameters  in  ordinary 
differential  equations  from  noisy  data  has  received  considerable 
attention  in  recent  years  as  indicated  by  the  large  number  of 
articles  [22,23,24,25,26],  review  papers  [27,28,29]  and  books  [10, 
30,31,32,33,34,35]  concerned  with  this  topic. 

Parameter  estimation  techniques  can  be  conveniently  classi¬ 
fied  as  being  either  sequential  or  nonsequential  methods.  The 
characteristic  feature  of  nonsequential  estimation  is  that  a  string 
of  input/output  data  is  collected  before  the  parameter  estimation 
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begins.  By  contrast,  in  sequential  estimation  a  parameter  estimate 
can  be  generated  at  each  sampling  instant.  Consequently,  sequential 
estimation  methods  are  well -suited  for  adaptive  control  applications 
where  process  parameters  change  and  are  to  be  "tracked"  (i.e.  esti¬ 
mated  on-line).  In  addition,  sequential  estimation  can  provide  esti¬ 
mates  of  state  variables  which  are  needed  for  control  but  often 
cannot  be  measured.  In  this  thesis,  state  variables  which  are 
impractical  or  inconvenient  to  measure  will  be  referred  to  as 
"inaccessible  state  variables". 

Sequential  estimation  of  states  and  parameters  is  often 
referred  to  as  a  filtering  problem  [29].  For  nonlinear  estimation, 
there  exist  several  nonlinear  filter  techniques,  such  as  the  extended 
Kalman  filter  [32,34],  Detchmendy  and  Sridhar  filter  [22],  Cox  filter 
[25,32]  and  many  minimal  variance  filters  [25,32].  These  methods 
can  be  used  for  filtering  out  noise  contained  in  the  process  and 
measurements,  and  estimating  unknown  parameters  and  inaccessible 
state  variables.  The  differences  between  these  nonlinear  filter 
techniques  lie  in  the  different  estimation  criteria  that  are  used 
and  the  approximations  that  are  involved  [35,19,30,31,32].  Since 
no  general  solution  of  the  filtering  problem  is  available,  all  of 
the  nonlinear  filter  techniques  require  approximations  to  be  made 
which  introduce  some  degree  of  suboptimality. 

Lately,  nonlinear  filter  schemes  have  been  considered  in 
chemical  engineering  problems.  Wells  [36]  used  the  extended  Kalman 
filter  to  perform  state  and  parameter  estimation  for  a  simulated 
single  reaction  taking  place  in  a  continuous  stirred  tank  reactor 
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(CSTR).  Gavalas  and  Seinfeld  [37]  used  a  Cox  filter  to  estimate 
catalyst  decay  rates  for  a  single  irreversible  reaction  taking  place 
isothermally  in  a  plug  flow  reactor.  The  extended  Kalman  filter  has 
been  used  by  Seinfeld  [16]  to  investigate  the  stochastic  features 
of  a  single  CSTR  reaction.  A  single  CSTR  control  scheme  with  time 
delay  in  the  feedback  control  loop  has  also  been  studied  [38]  using  a 
Detchmendy  and  Sridhar  nonlinear  filter  to  estimate  the  unknown 
kinetic  parameter  of  a  single  reaction.  In  a  batch  reactor  applica¬ 
tion,  Lee  [30]  used  a  Detchmendy  and  Sridhar  nonlinear  filter  to 
estimate  the  unknown  rate  constant  of  a  single  isothermal  reaction. 
Pell  [39]  used  a  different  nonlinear  filter  to  estimate  the  unknown 
kinetic  parameters  of  a  simple  consecutive  reaction  for  both  iso¬ 
thermal  and  nonisothermal  conditions.  Coggan  and  Noton  [40]  have 
used  the  extended  Kalman  filter  to  investigate  the  filter  performance 
by  simulating  a  complicated  blending  process  and  a  thermal  system. 
Applications  of  the  proposed  filtering  schemes  to  physical  systems 
of  interest  to  chemical  engineers  have  seldom  been  reported.  A  note¬ 
worthy  exception  is  the  application  of  the  extended  Kalman  filter  to 
an  industrial  reactor  system  [41].  However,  few  details  of  this 
application  are  available. 

With  so  many  nonlinear  filters  available  in  the  literature  and 
since  different  filters  have  been  applied  successfully  in  different 
examples,  comparisons  of  different  filters  have  also  been  made 
[23,25,32].  Wishner  et  al  [23]  compared  the  performance  of  the 
extended  Kalman  filter,  a  second  order  nonlinear  filter  and  a  single 
stage  iteration  filter,  which  have  the  same  theoretical  basis  but  are 
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successively  more  complex  to  implement.  From  their  simulated  results 
for  three  numerical  examples,  they  concluded  that  the  performance  of 
the  filter  is  successively  improved  by  the  additional  complexity. 
Schwartz  and  Stear  [25]  compared  the  performance  of  eight  different 
filters,  one  linear  and  seven  nonlinear,  on  two  scalar  nonlinear 
systems.  They  concluded  that  none  of  the  eight  filters  provided  a 
clearcut  superiority  over  the  others.  Jazwinski  [32]  later  criticized 
their  conclusions  on  the  ground  that  the  large  noise  levels  in  the 
Schwartz  and  Stear 's  examples  masked  the  nonlinearity  of  the  measure¬ 
ment  equation;  hence,  the  simulated  filters  in  their  examples  cannot 
be  recommended  on  the  basis  of  superior  performance.  However,  he 
then  concluded  that  the  relative  value  of  different  nonlinear  filters 
in  a  particular  problem  cannot  be  assessed  a  priori,  but  must  be 
determined  by  simulation. 

In  this  thesis,  the  Detchmendy  and  Sridhar  nonlinear  filter 
[22]  will  be  applied  to  two  batch  reactor  examples.  Since  their 
filter  is  based  on  a  weighted  least  square  performance  index,  no 
prior  statistical  information  about  the  process  noise  is  required. 
However,  the  weighting  matrices  used  in  the  performance  index  must  be 
speci fied. 

[221 

2.2  Summary  of  Detchmendy  and  Sridhar  Filter1-  J  Algorithm 

In  this  section,  Detchmendy  and  Sridhar's  filter  algorithm  [22] 
is  briefly  described.  A  more  detailed  derivation  is  presented  in 
Appendix  A. 

Consider  a  nonlinear  system  governed  by  the  ordinary  vector 
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differential  equation 

z(t)  =  £[z_(t),£,t]  (2-1) 


where  _z  and  £  are  vectors  of  state  variables  and  time  invariant  para¬ 
meters,  respectively,  and  t  denotes  time.  Define 


x  = 


then  Equation  (2-1)  can  be  rewritten  as 


X  = 


1 

0 


or 


x  =  f[x(t),t] 


(2-2) 


(2-3) 


(2-4) 


Let  the  total  dimension  of  z  and  £  be  n,  then  x_  becomes  a  n-dimensional 
vector  and  _f (x_,t)  is  a  n-dimension  vector  function. 

In  general  for  physical  systems,  dynamic  errors  (e.g.  random 
noise,  modelling  errors)  are  often  involved  in  the  process.  Let  the 
dynamic  error  present  in  the  process  be  represented  by  an  additive 
term  in  the  vector  differential  equation.  Then  the  dynamic  model  of 
the  process  can  be  expressed  by 

x(t)  =  f(x,t)  +  K(x.,t)  w(t)  (2-5) 

and  let  the  observation  model  also  contain  additive  noise,  ^(t) ,  with 
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X(t)  =  h_(x,t)  +  v(t)  (2-6) 

where 

y_(t)  =  m-vector  of  outputs  (i.e.  measured  variables) 

K(x,t)  =  n  xi  matrix,  h,(x,t)  =  m-vector  function 

w(t)  =  ^-vector  of  unknown  y_{t)  =  m-vector  of  measurement 

inputs  ,  errors 

with  w(t)  and  y_(t)  assumed  to  be  random  variables  with  unknown 
statistics . 

By  measuring  y_(t),  it  is  desired  to  find  the  least  squares  esti¬ 
mate  of  _x(t),  denoted  by  x(t),  such  that  the  following  performance 
index  is  minimized. 


I  = 


C||y.(t)  -  h(x,t)|||(t)  +  II x.  -  f(i.t)||^tj]dt  (2-7) 


where  t^  is  the  specified  final  time  and  Q_(t)  and  W(t)  are  weighting 
matrices  which  determine  the  relative  weighting  to  be  placed  on  the 
individual  terms  in  the  performance  index. 

By  applying  Pontryaginls  maximum  principle  to  the  dynamic 
optimization  problem  posed  by  Equations  (2-4) - (2-7) ,  and  making 
certain  approximations (see  Appendix  A),  Detchmendy  and  Sridhar  [22] 
derived  the  following  filter 

ajl' (x,t) 

x(t)  =  f(x,t)  +  £.(  t )  [ - - - ]Q.(t)[y_(t)-h(x,t)]  (2-8) 

-  9x  “ 


I 
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8f(x,t)  3  f 1 (x , t ) 

P(t)  =  [ - 7 - ]P(t)  +  P(t)[ - =—] 

~  dX  ~  ~  3X 

3h 1 (x  ,t) 

+  p(t)  —  [(— -1  )a(t)(x(t)  -  k(x,t))]  p(t) 

“  3X  3X.  “  ~ 

+  K(x,t)  V_1(x,t)  K'(x,t)  (2-9) 

with 


V(x,t)  =  £(x,t)  K(t)  £(x,t)  (2-10) 

where  h'  (x_,t)  and  K_'(x_,t)  are  the  corresponding  transposes  of 
h_Cx,t)  and  K(x_,t) ,  V_”^  (x_,t)  is  the  inverse  of  V_(x_,t)  and  P_(t )  is  an 
nxn  symmetric  matrix  which  can  be  roughly  interpreted  as  the  co- 
variance  of  the  random  vector  x_(t)  [37],  Integration  of  Equations 
(2-8)  and  (2-9)  provides  an  estimate  _x(t)  of  the  actual  state  .x(t) 
based  on  measurements  of  the  output,  y_(t)  - 

In  this  filter  estimation,  _x  and  are  always  coupled  and 
evaluated  simultaneously  after  specifying  the  initial  values  of 
Q(0),  W(0),  P_(0)  and  _x(0)  and  choosing  a  final  time  t^„  Hence  the 
total  number  of  differential  equations  that  are  investigated  during 
filtering  becomes  n  +  n(n  +  l)/2  =  n(n  +  3)/2„  By  analogy  with  the 
Kalman  filter  which  is  optimal  for  a  linear  system,  the  weighting 
matrices,  Q_(t)  and  W(t),  can  be  chosen  as  the  inverse  of  the  co- 
variance  matrices  of  the  output  measurement  error,  v_(t) ,  and  input 
noise,  K(x,t)w(t),  respectively  [31].  Similarly,  for  the  Kalman 
filter,  P_(0)  is  chosen  to  be  the  inverse  of  the  covariance  matrix 
for  the  initial  state,  x(0),  assuming  it  is  a  random  variable  with 
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unknown  statistics  [31,37].  Since  approximations  are  used  in  deriv¬ 
ing  the  filter  equations,  the  initial  value  of  p_(0)  is  not  known 
exactly;  hence  numerical  experimentation  is  usually  required  to  find 
a  satisfactory  value  [31], 

2.3  Simplified  Filter  Algorithm  -  I 

In  using  Detchmendy  and  Sridhar's  nonlinear  filter,  -n-4|— 
differential  equations  must  be  integrated  to  determine  n  state 
variables  (including  the  unknown  parameters).  Since  the  estimation 

A 

of  _x  and  P_  are  always  coupled  and  solved  simultaneously  at  each 
sampling  instance,  the  filter  estimation,  which  requires  matrix 
manipulation  and  integration,  becomes  very  time  consuming  when  n  is 
large.  However,  in  most  industrial  situations,  the  digital  computer 
used  for  real-time  applications  is  a  small  or  medium  size  computer, 
such  as  the  IBM  1800,  with  restrictions  on  memory  size  and  execution 
speed.  Thus  the  computation  time  of  filter  estimation  becomes  a 
severe  point  for  real-time  applications. 

For  the  extended  Kalman  filter,  Dressier  and  Ross  [42]  have 
stated  that  they  investigated  several  approximate  filter  algorithms 
in  order  to  reduce  the  computation  time  without  significantly 
degrading  the  filter  performance.  They  found  the  most  promising 
simplified  filter  algorithm  to  be  one  which  uses  a  piecewise 
recursive  estimation  [42],  The  basic  idea  of  the  proposed  piecewise 
recursive  estimation  scheme  is  to  separate  the  filter  estimation  into 
two  estimation  "loops";  one  for  state  estimation  and  one  for  the 
covariance  equations.  This  can  be  explained  as  follows. 


- 
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In  a  typical  filtering  application,  the  state  estimate  x(t) 
and  the  covariance  matrix  P_(t)  are  generated  by  numerical  integration 
using  the  data  sampling  interval,  t  ,  as  the  integration  interval. 

Thus  a  sequence  of  values,  x(nt  )  and  P(n t  ),  are  generated  for 
n  =  0,1,2,...  Dressier  and  Ross  [42]  proposed  that  the  covariance 
equations  be  updated  every  M  sampling  intervals  where  M  is  an 
integer  greater  than  one.  This  strategy  will  result  in  a  significant 
reduction  in  the  computational  requirements  for  the  filter.  It  also 
has  the  effect  of  maintaining  P_(t)  constant  over  M  sampling  interval 

/V 

while  _x  is  still  updated  every  sampling  interval.  In  effect  P£t) 
becomes  piecewise  constant  and  the  amount  of  computations  required 
to  integrate  the  covariance  equations  is  reduced  by  a  factor  of  M. 

Dressier  and  Ross  [42]  have  applied  their  "piecewise  recursive" 
modification  of  the  extended  Kalman  filter  to  a  simple  second  order 
numerical  example.  Apparently,  this  is  the  only  published  applica¬ 
tion  of  this  approach.  In  this  study,  the  piecewise  recursive  modi¬ 
fication  will  be  further  evaluated  as  a  modification  to  the  Detchmendy 
and  Sridhar  filter  and  by  application  of  the  resulting  filter  to  two 
batch  reactor  examples. 

2.4  Simplfied  Filter  Algorithm  -  II 

A  second  type  of  simplified  nonlinear  filter  technique  has 
been  proposed  by  Seinfeld  and  Chen  [14].  The  basic  idea  is  to  replace 
the  covariance  differential  equations  by  algebraic  expressions  for 
the  elements  of  P_(t).  Since  the  computation  time  required  to  inte¬ 
grate  differential  equations  will,  in  general,  be  much  greater  than 
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that  for  evaluating  simple  algebraic  expressions,  the  reduction  in 
computational  time  and  storage  requirements  can  be  very  significant 
for  real  time  appl ications .  A  systematic  method  for  determining  the 
algebraic  expressions  is  to  examine  the  time  behavior  of  _P(t)  during 
estimation  using  the  original  filter.  Then  the  response  of  the 
elements  of  P(t)  can  be  roughly  approximated  by  the  use  of  simple 
functions  such  as  exponentials  [14]. 

The  only  known  application  of  this  technique  is  the  reactor 
example  in  the  unpublished  work  by  Seinfeld  and  Chen  [14].  Hence 
further  evaluation  of  this  approach  would  appear  to  be  desirable. 


17 


CHAPTER  III 

STUDIES  OF  A  SINGLE  REACTION 


3.1  Formulation 

The  optimal  temperature  control  for  a  chemical  reaction  occur¬ 
ring  in  a  batch  reactor  has  two  equivalent  senses:  it  will  achieve 
the  maximum  yield  in  a  fixed  time  and  it  will  require  the  minimum 
amount  of  time  to  achieve  a  certain  conversion.  For  a  single  reaction, 
only  reversible  exothermic  reactions  provide  interesting  optimization 
and  control  problems.  By  contrast,  in  endothermic  reactions  heat  is 
absorbed  and  the  optimal  control  policy  is  to  provide  as  much  heat  as 
is  possible.  For  a  single  exothermic  reaction  the  optimal  temperature 
control  policy  is  to  maximize  the  reaction  rate  at  each  instant  of 
time  [3,11].  It  has  been  found  that  the  optimal  temperature  of  a 
single  reversible  exothermic  reaction  is  somewhat  less  than  the 
equilibrium  temperature  at  any  given  composition  [3,11]. 

Consider  a  single  reversible  exothermic  reaction  taking  place 
in  a  constant  volume  batch  reactor.  If  the  chemical  species  present 
are  A-j  ,  A2»  ...  A$  the  reaction  stoichiometry  can  be  written  as  [3] 

l  a.  A.  =0  (3-1) 

1-1  1  1 

where  a.,  the  stoichiometric  coefficient  of  A., is  chosen  to  be  posi¬ 
tive  for  a  product  and  negative  for  a  reactant. 

If  C.  is  the  concentration  of  A.,  at  any  time  t,  and  CiQ  its 
concentration  at  time  t  =  0,  then  the  extent  of  reaction  c  is  defined 


' 
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as  [3] 


c  E  (Cio  “  Ci)/ai  (3-2) 

During  the  reaction,  two  governing  equations,  namely  the  mass 
and  energy  balances,  can  be  used  to  describe  the  state  of  the  batch 
reactor.  These  well-known  equations  can  be  written  as  [3,11] 

g§  =  r(c,T)  (3-3) 

Cp  31  =  (-AH)r(c.T)  -  q'  (3-4) 

where  r(c,T),  the  reaction  rate,  is  a  function  of  c  and  temperature  T, 
C  is  the  total  heat  capacity  per  unit  volume  (including  inerts), 

r 

(aH)  is  the  heat  of  reaction  (negative  for  an  exothermic  reaction), 
and  q 1  is  the  rate  of  heat  removal. 

If  it  is  assumed  that  C  is  constant,  then  Equation  (3-4)  can 

r 

be  written  as : 

gg  =  J  r(c,T)  -  q  (3-5) 

where  J  =  (-AH)/C  and  q  =  q'/Cp 

A  common  form  of  the  reaction  rate  expression  [3]  is 

s  8.  s  y, 

r(c,T)  =  k-!exp(-E-,/RT )  n  c.  -  kAexp(-E?/RT)  n  c,  (3-6) 

l  1  i=l  1  L  L  i=l 

where  k!  is  the  frequency  factor,  E.  is  the  activation  energy,  R  is 
the  gas  constant  and  6-  and  are  the  forward  and  backward  reaction 


orders . 
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Since  in  the  optimal  control  policy,  the  instantaneous  reaction 
rate  is  a  maximum,  the  optimal  temperature  schedule  will  be  the  solu¬ 
tion  of 


ar(c,T)  . 
3T 


0 


(3-7) 


Thus  differentiating  Equation  (3-6)  with  respect  to  T,  gives 
the  following  optimal  temperature  sequence  as  a  function  of  c  [3] 


T  (c) 
m 


E2  -  E] 


RTln(k^r')  n  c.s 

1  =  1  1 


(3-8) 


where  £..  =  y..  -  3^  and  T  (c)  denotes  the  optimal  temperature  sequence. 

Substituting  Equation  (3-8)  into  Equation  (3-6),  the  correspond¬ 
ing  maximum  reaction  rate  is 


rm[c’Tm(c)]  = 


S  3i  P+1 

{[ki/(p+l)] .n  c.  } 

_  I _ 1  =1  I _ 

s  Yi  p 
[Ck'/p)  .n  c,  ’] 

1  =  1  1 


(3-9) 


where 


p  =  E-,/(E2  -  E,) 


(3-10) 


Then  from  Equation  (3-4)  the  rate  of  heat  removal  at  extent  c 
to  keep  the  reactor  on  the  optimal  path  is 


■  cp  “  U-AH)r[c,Tm(c)]}-  Cp  ^ 


(3-11) 


But  on  the  optimal  path 
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so  that 


dT  =  dVc)  dc 
dt  dc  dt 


T'(c)  r  [c,T  (c)] 
m  m  m v  J 


0-12) 


■  Cp[J-r;(c)]  r^.yc)]  (3-13) 

Equations  (3-8),  (3-9),  and  (3-13)  will  provide  a  straightforward 
computation  for  the  optimal  temperature  schedule  and  the  optimal  cool¬ 
ing  rate  as  a  function  of  c. 

\ 

3 . 2  Chien  and  Aris  Example 
3.2.1  Optimal  Control  Policy 

Consider  the  hypothetical  second  order  reversible  reaction 
studied  by  Chien  and  Aris  [5,6]  with  reaction  mechanism 

A]  +  A2  t  A3  (3-14) 

and  reaction  rate  given  by 

r(c,T)  =  k-JexpC-E-j/RT)  (3-c)  (2-c)  -  k^exp(-E2/RT)  (1+c)2  (3-15) 

The  kinetic  parameters  used  by  Chien  and  Aris  [5,6]  are: 

k'  =  exp (12)  ,  k2  =  exp (30) 

E]  =  24000  ,  ( — AH )  =  E2-E1  =  24000 

J  =  (-AH)/C  =  400  p  =  1 

r 

From  Equations  (3-8),  (3-9)  and  (3-13),  the  optimal  temperature 


schedule  is, 


. 


the  maximum  reaction  rate  is  given  by, 


[(k'/2)(3-c)(2-c)]2 
40  +  c)2 


(3-17) 


and  the  optimal  cooling  rate, 


(3-18) 


where 


(3-19) 


Although  the  optimal  temperature  and  cooling  rate  computations 
are  quite  straightforward,  it  should  be  noted  that  for  small  values 
of  c  the  optimal  temperature  will  be  very  high  and  may  exceed  the 
limits  of  safety  or  may  require  an  excessively  large  cooling  rate. 

When  the  initial  state  of  the  reactor  is  known  and  only  cooling  (but 
not  heating)  capacity  is  available,  it  has  been  shown  [4,9]  that  the 
optimal  control  policy  is  a  combination  of  an  optimal  cooling  path 
and  either  an  adiabatic  path  or  a  path  with  maximum  cooling,  depending 
upon  the  initial  state  of  the  reactor  [9],  The  optimal  control  policy 
is  shown  in  Figure  3-1  where  the  reactor  is  initially  operated 
adiabatically  and  then  switched  to  the  optimal  cooling  path  after  a- 


. 
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certain  temperature  is  attained  [5,43].  This  can  be  explained  as 
follows . 

For  adiabatic  operation,  q1  =  0  and  Equation  (3-5)  becomes 

^  =  J  r(c,T)  (3-20) 

From  Equations  (3-3)  and  (3-20)  it  follows  that 

c  =  c0  +  (T  -  T0)/J  (3-21) 

where  T  is  the  initial  temperature  when  the  extent  is  c  . 

o  o 

If  along  the  adiabatic  path,  T  is  measured  and  the  expression 
in  Equation  (3-22)  is  evaluated, 

Tm[co  +  (T  ‘  To)/Jl  '  T  (3-22) 

then  the  switch  to  the  optimal  cooling  path  should  be  made  as  soon  as 
Equation  (3-22)  becomes  zero  (i.e.  as  soon  as  the  reaction  path 
intersects  the  switching  curve). 

In  Equation  (3-18),  the  optimal  cooling  rate  is  a  function  of 
c.  However,  in  many  reactions  the  extent  cannot  be  directly  measured 
while  the  reaction  is  taking  place.  Since  temperature  measurement  is 
easier,  cooling  rate  based  on  temperature  (rather  than  extent)  will  be 
more  convenient.  It  is  known  [6]  that  along  the  switching  curve,  the 
optimal  temperature  is  monotone  decreasing;  consequently.  Equation 
(3-16)  can  be  inverted  and  c  obtained  as  a  function  of  T.  Substitu¬ 
tion  of  the  resulting  expression  for  c  into  Equation  (3-18)  gives  the 
optimal  cooling  rate  q^  as  a  function  of  temperature. 
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An  interesting  feature  which  characterizes  the  optimal  control 
path  is  the  inherent  self-stability  when  temperature  measurement  is 
used  [6].  Suppose  the  reaction  state  is  off  the  optimal  path  at  point 
B  in  Figure  3-2  where  the  extent  is  greater  than  that  at  the  optimal 
point,  A.  If  the  optimal  rate  of  cooling  based  on  the  measured  tempera¬ 
ture  is  applied,  then  the  state  tends  to  change  toward  the  optimal 
path.  This  occurs  since  the  reaction  rate  at  B  is  less  than  that  at 
the  corresponding  optimal  point,  A,  so  that  the  rate  of  cooling  is 
greater  than  the  rate  of  heat  generation.  On  the  other  hand,  if  the 
reaction  state  is  at  point  C  where  the  extent  is  less  than  the 
corresponding  optimal  point  A  and  the  optimal  cooling  rate  based  on 
point  A  is  applied,  then  the  heat  generation  will  preponderate  and 
again  will  drive  the  reactor  state  to  the  optimal  path.  By  contrast 
if  measurements  of  c  are  used  to  calculate  the  optimal  cooling  rate, 
the  effect  will  be  to  drive  the  trajectory  away  from  the  optimal 
path  [6]. 
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Figure  3-1  Optimal  Operation  for  a  Single  Reaction 
Taking  Place  in  a  Batch  Reactor  [6] 


Figure  3-2  Stability  of  the  Optimal  Path  [6] 
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3.2.2  Parameter  Sensitivity 

It  is  well-known  that  the  optimal  control  policy  is  based  on 
the  availability  of  a  precise  mathematical  model.  But  in  real 
physical  situations,  the  formulated  mathematical  model  is  only 
approximate  and  the  process  parameters  are  difficult  to  determine, 
particularly  for  chemical  reaction  systems.  To  see  how  inaccurate 
parameter  values  affect  the  system  performance,  a  parameter  sensi¬ 
tivity  analysis  will  be  made  to  assign  accuracy  requirements  for  a 
system  model  and  to  give  a  deeper  insight  into  control  system  design. 

In  the  previous  discussion,  the  optimal  control  policy  was 
implemented  by  measuring  temperature  along  the  adiabatic  path  and 
switching  to  the  optimal  cooling  rate  after  Equation  (3-22)  became 
zero  or  attained  a  preset  small  value.  To  determine  the  effect  of 
activation  energy  accuracy  on  the  control  policy,  errors  of  -  10% 
in  activation  energy  E-j  were  used  in  calculating  the  optimal  cool¬ 
ing  rate.  In  this  simulation,  the  initial  reaction  conditions  are 

A 

c(0)  =  0  and  T (0 )  =  500°K,  and  are  assumed  to  be  known.  The  result¬ 
ing  reactor  trajectories  are  shown  in  Figure  3-3  when  either 
temperature  measurements  or  extent  measurements  are  used  to  calculate 
the  optimal  cooling  rate.  The  effect  of  using  inaccurate  values  of 
E-j  (and  temperature  measurements)  to  calculate  the  optimal  cooling 
rate  is  illustrated  in  Table  3-1.  For  a  +  10%  error  in  the  assumed 
value  of  E-,,  the  time  to  completion  is  8050  sec.  which  is  much  greater 
than  the  optimal  value  of  3880  sec.  Consequently,  the  potential  value 
of  on-line  parameter  estimation  in  reducing  the  time  to  completion  is 
readily  apparent. 
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Figure  3-3  Error  in  Optimal  Path  Introduced  by  Calculating  the  Cooling  Rate  with  -  10%  Error 
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TABLE  3-1 

EFFECT  OF  ERROR  IN  THE  ASSUMED  VALUE  OF 

E-j  ON  TIME  TO  COMPLETION  (WHERE  c  =  1.4) 


Error  in  E-j  Time  to  Completion  (sec) 


0 

3,880 

+  1% 

3,913 

-  1% 

3,923 

+10% 

8,050 

-10% 

16,360 

*  Based  on  nominal  value  of  E-j  =  24,000 

3.2.3  Suboptimal  Control  Policy 

One  of  the  limitations  of  most  optimal  control  formulations  is 
the  assumption  that  the  system  parameters  are  known  and  time  invariant. 
In  many  practical  situations,  this  is  not  true.  Even  when  nominal 
parameter  values  are  known,  they  may  be  subject  to  change.  For 
example,  in  batch  reactions,  the  catalyst  activity  may  change  during 
the  reaction  or  undetected  changes  in  initial  charge  characteristics 
may  occur  from  batch  to  batch.  Consequently,  when  the  optimal  control 
policy  is  quite  sensitive  to  parameter  variations,  control  based  on 
inaccurate  parameter  values  will  give  poor  process  performance. 

In  chemical  engineering  applications,  several  authors  have 
proposed  feedback  control  schemes  which  are  suboptimal  and  attempt 
to  meet  two  main  objectives  [13,14,15,16].  One  of  the  objectives  is 
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to  avoid  the  complicated  computations  required  to  calculate  the  optimal 
control  policy,  the  other  is  to  compensate  for  parameter  changes  and 
process  environmental  variations. 

Here  the  conventional  proportional  control  scheme  proposed  by 
Mi  liman  and  Katz  [13]  will  be  investigated  for  the  Chien  and  Aris 
reaction  example.  With  temperature  as  the  manipulated  variable,  the 
proportional  control  policy  [13]  is  given  by 

T  =  K-,  +  K2  c  (3-23) 

where  K-j  and  are  constants. 

By  using  a  gradient  method  and  repeated  integration  of  the  mass 
balances,  Mi  liman  [12^13]  determined  the  best  values  of  K-j  and  by 
assuming  the  kinetic  parameters  were  known  exactly  and  did  not  vary 
with  time.  His  results  indicate  that  for  several  batch  reaction 
examples,  the  proportional  control  scheme  resulted  in  a  time  to 
completion  which  was  very  close  to  the  optimal  value  .  However,  the 
sensitivity  of  the  proportional  control  scheme  to  variations  in  the 
assumed  kinetic  parameters  was  not  investigated. 

In  the  Chien  and  Aris  example,  the  nonadiabatic  temperature 
schedule,  Tm(c),  is  almost  a  straight  line  (see  Figure  3-3).  Thus, 
approximate  values  of  K-j  and  can  be  determined  graphically  to 
give  the  following  suboptimal  control  law 

T  =  702  -  84  c  (3-24) 

Simulation  results  showed  that  if  Millman's  control  scheme  in 
Equation  (3-24)  was  initiated  at  the  adiabatic  switching  point,  a 
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completion  time  of  3885  sec.  would  result.  This  value  is  very  close 
to  the  optimal  value  of  3880  sec.  as  indicated  in  Table  3-1.  If  the 
extent  of  reaction  is  measurable  or  can  be  estimated,  this  control 
policy  seems  to  be  quite  effective. 

Optimal  cooling  rate  based  on  an  incorrect  value  of  the 
activation  energy  E-j  can  result  in  poor  reactor  control  as  has  been 
shown  in  Section  3.2.2.  Hence,  cooling  rate  based  on  an  optimal 
policy  is  not  feasible  when  changes  in  activation  energy  E-j  occur 
and  are  not  detected.  To  determine  how  the  proportional  control 
policy  works  when  the  actual  (rather  than  the  assumed)  value  of  E-j 
changes,  simulated  results  of  optimal  and  proportional  control 
strategies  are  shown  in  Table  3-2.  It  is  remarkable  that  the  pro¬ 
portional  control  policy  in  Equation  (3-24)  results  in  completion 
times  which  are  within  1%  of  the  optimal  values  despite  the 
occurrence  of  10%  changes  in  E^ .  By  contrast,  the  sensitivity  of 
the  optimal  control  policy  to  parameter  variation  is  again  apparent 
from  Table  3-2.  In  these  optimal  control  calculations  (and  the 
ones  to  follow),  the  cooling  rate  is  calculated  from  temperature 


measurements . 
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TABLE  3-2 

EFFECT  OF  VARIATIONS  IN  THE  ACTUAL  VALUE  OF 

E-j  ON  CONTROL  SCHEMES 

Control  Variation  in  E-j  Time  to  Completion(sec) 


optimal  based  on  actual 
value  of  E-j 

-10% 

473 

11 

+10% 

32,685 

optimal  based  on  nominal 
value  of  E-j 

-10% 

1  ,046 

II 

+10% 

132,123 

proportional  based  on 
nominal  value  of  E-j 

-10% 

473 

II 

+10% 

32,830 

*  Based  on  the  nominal  value  of  E-j  =  24,000 

Thus  it  appears  that  the  proportional  control  policy  can  com¬ 
pensate  for  variation  of  parameters  even  when  the  changes  are  not 
detected.  However,  if  the  extent  c  cannot  be  readily  measured,  some 
type  of  on-line  state  estimation  would  be  required.  This  topic  is 
considered  in  the  next  section, 

3.3  Sequential  State  and  Parameter  Estimation 
3.3.1  Parameter  Estimation 

Optimal  and  proportional  control  policies  have  been  investi¬ 
gated  for  the  Chien  and  Aris  reaction  example.  Parameter  sensitivity 
calculations  showed  that  an  incorrect  value  of  E-j  will  result  in  poor 


. 
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reactor  control  if  the  optimal  cooling  policy  is  used.  In  Millman's 
method  [12,13],  the  calculation  of  constants  K-j  and  in  Equation 
(3-23)  is  based  on  assumed  process  parameter  values,  and  is  quite 
insensitive  to  changes  in  E-j  „  However,  Millman's  method  requires 
that  an  estimate  (or  measurement)  of  the  extent  of  reaction  must  be 
avai 1  able . 

With  reactor  operation  divided  into  an  adiabatic  path  and  a 
optimal  path  as  shown  in  Figure  3-1,  Chien  and  Aris  [5]  used 
numerical  approximation  of  the  time  derivative  of  temperature  (i.e. 
dT/dt  in  Equation  (3-20))  and  matrix  inversion  least  squares  tech¬ 
nique  to  identify  the  unknown  parameters  during  the  adiabatic  path. 

Ray  and  Aris  [7]  used  the  same  technique  to  identify  unknown 
parameters  from  several  isothermal  runs  then  switched  to  a  non- 
isothermal  path  (i.e.  optimal  path  in  Figure  3-1),  This  latter 
approach  will  result  in  some  degree  of  suboptimality  in  the  whole 
process.  One  common  problem  to  both  approaches  is  the  large  round¬ 
off  computation  errors  encountered  as  the  determinant  of  the  matrix 
to  be  inverted  becomes  very  small.  Their  methods  have  the  advantage 
that  initial  guesses  of  the  parameters  are  not  required.  However, 
their  methods  also  suffered  two  main  shortcomings.  One  is  that  the 
computation  scheme  cannot  easily  be  extended  to  a  complex  reaction 
scheme,  since  in  this  case  nonlinear  regression  techniques  would 
normally  be  required.  The  other  disadvantage  is  that  the  accuracy 
of  the  required  data  is  quite  high  so  when  the  measurement  noise  is 
high  (even  for  a  measurement  error  of  1°C),  their  methods  will  not 
provide  satisfactory  estimates. 
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As  an  alternative  approach,  the  Detchmendy  and  Sridhar  nonlinear 
filter  technique  [22]  can  be  used  to  perform  parameter  estimation 
during  the  adiabatic  path.  Suppose  E-j  is  the  only  undetermined  para¬ 
meter  in  this  example  and  assume  no  dynamic  errors  are  involved  in  the 
process.  It  will  also  be  assumed  that  temperature,  but  not  extent, 
can  be  measured  and  that  the  noisy  temperature  data,  y(t),  can  be 
represented  as 


y(t)  =  T(t)  +  e (0  ,a ) 


(3-25) 


where  the  measurement  errors  e(0,a)  are  assumed  to  be  Gaussian 
distributed  with  zero  mean  and  an  appropriate  standard  deviation,  a. 

One  of  the  difficulties  involved  in  parameter  estimation  by  a 
nonlinear  filter  technique  is  that  the  filter  differential  equations 
are  often  stiff,  i.e.  the  eigenvalues  differ  greatly  in  magni tude[4£, 
14].  To  reduce  this  difficulty,  Seinfeld  and  Chen  [14]  proposed 
defining  new  variables,  containing  the  unknown  parameters,  which  are 
of  the  same  order  of  magnitude  as  the  states.  To  reduce  the  stiff¬ 
ness  effect  in  this  example,  the  state  variables  of  the  filter  are 
defined  as 


x,(t) 

x2(t) 

x3(t) 


c(t) 

T ( t)  -  500 

100 

M1)  ~  E1  ,ref 

E1  ,ref 


(3-26) 

(3-27) 

(3-28) 


where  c(t)  is  the  estimated  extent,  E] (t)  and  f(t)  are  the  estimated 
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activation  energy  and  temperature,  and  E-j  ^  is  a  specified  refer¬ 
ence  value. 

Then  the  filter  equations  (2-8  to  2-10)  can  be  developed  with 
the  following  state  functions  and  Jacobian  applicable  during  the 
adiabatic  path: 


q  =  k|exp[-(E-|  (t)/R)/T (t)](3-x-| )  (2-^ ) 

E,(t)+24000  ,  , 

-  k£exp[-(— - p - )/T(t)] (1+X-, (3-29) 


f2  =  400  f] 


f3  ‘  0 


(3-30) 

(3-31) 


=  k-jexp[-(E-|  Ct)/R)/T (t) ] (-5  +  2x-j ) 


E, (t)+24000 

-  2k2exp[-(— - ^ - )/T (t)](l+x-j ) 

A. 

f12  =  { k-j exp[- (E-j  (t)/R)/T(t)](— )(3-x-|)(2-x-| ) 


(3-32) 


E,(t)+24000  E,  (t)+24000  ,  ?  lnn 

-  k,exp[.(^__ - -)/fCt)]t-! - r - ■)(l«1)ZlA)  (3- 


f13  =  {qexp[-(E1(t)/R)/T(t)](3-x1)(2-x1) 


E,(t)+24000  „  „  o  E,  f 

k~exp[-(— - 5 - )/T(t)](l+x1)2H-  -rf^) 

2  R  1  RT(t) 


(3-34) 


33) 
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f21 

400 

fll 

(3-35) 

f22  = 

400 

fl  2 

(3-36) 

f23  = 

400 

f  1 3 

(3-37) 

f31  = 

f32 

=  f33  =  0 

(3-38) 

with 

^(t) 

=  E1 

,ref  +  x3(t)  E1  ,ref 

(3-39) 

f(t)  =  500  +  100  x2(t)  (3-40) 

Suppose  the  initial  reactor  state,  c  =  0  and  T  =  500  at  t  =  0, 
is  known  exactly  and  E]  ref  is  chosen  to  be  24000  (the  actual  value 
of  E-j ) ;  then  if  the  initial  estimate  of  E-,  is  10%  low  (i.e. 

E-j(O)  =  21600),  the  initial  states  of  the  filter  are 

x]  (0)  =  0  ,  x2(0)  =  0,  x3(0)  =  -0.1 

As  discussed  in  Section  2.2,  numerical  experimentation  is 
required  to  determine  a  satisfactory  initial  covariance  matrix  £(0). 
For  the  first  trial  run,  the  matrix  P_(0 )  was  arbitrary  set  equal  to 
the  identity  matrix.  Additional  trials  indicated  that  a  diagonal 
P_(0)  matrix  was  adequate  providing  that  the  element  corresponding  to 
the  unknown  parameter  is  chosen  to  be  larger  than  the  corresponding 
diagonal  elements  for  the  states.  After  10  trials,  the  following 
initial  covariance  matrix  £(0)  was  finally  selected 
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0.0025 


0 

0.01 

0 


0 

0 


P(0) 


0 

0 


Thus  with  P_(0)  specified  and  the  initial  state  available, 
simulation  of  the  process  estimation  was  performed  by  integrating 
the  filter  equations  with  a  Runge-Kutta-Gi 1 1  4th  order  routine  [43] 
using  a  data  sampling  interval  of  t  =  1 .  The  process  data  were 
generated  by  integrating  Equations  (3-5)  and  (3-25)  with  the  same 
integration  routine  and  an  integration  step  size  of  0.5  sec.  The 
noise  data  was  generated  by  a  standard  Gaussian  random  number  genera¬ 
tor  with  a  standard  deviation  of  as  0.1  and  zero  mean. 

Figure  3-4  shows  the  estimate  of  E-j  obtained  using  the  standard 

A 

full  filter  (Filter  I)  when  the  assumed  values  of  E-j(O)  contain  10% 
errors.  The  parameter  estimation  is  judged  to  be  quite  satisfactory 
since  convergence  was  obtained  before  the  switching  time  of  1695 
sec.;  i.e.  after  t  >  1695  sec.  the  adiabatic  operation  is  replaced  by 
optimal  cooling  which  requires  an  estimate  of  E-j . 

In  order  to  evaluate  the  performance  of  the  piecewise  recur¬ 
sive  filter  algorithm  (Filter  II)  as  discussed  in  Section  2.3,  the 
update  interval  for  the  covariance  differential  equation  was  chosen 
to  be  tp  =  10  (i.e.  10  times  the  data  sampling  interval).  Then  with 
the  same  initial  state  and  initial  covariance  matrix  P_(0)  as  used  in 
Filter  I,  the  process  is  simulated  under  the  same  conditions  of 
Figure  3-4.  Results  are  shown  in  Figure  3-5.  In  comparison  with 
Figure  3-4  for  the  Filter  I  estimation,  the  Filter  II  estimation 
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Figure  3-5  Filter  II  Estimation  of  E-,  with  -  10%  Error  in  E-,(0)(t  =  1  sec. 
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shows  a  more  rapid  convergence.  Hence  the  Filter  II  algorithm  appears 
to  be  a  promising  estimation  scheme. 

As  discussed  in  Section  2=4,  the  computation  time  required  for 
the  nonlinear  filter  technique  is  a  critical  factor  in  real-time 
applications  because  integration  of  a  large  number  of  differential 
equations  is  very  time-consuming.  If  the  covariance  differential 
equations  are  replaced  by  algebraic  equations,  then  the  saving  in 
computation  time  can  be  quite  significant.  In  choosing  appropriate 
algebraic  equations,  the  time  variation  of  each  element  of  the  co- 
variance  matrix  using  Filter  I  or  Filter  II  was  first  obtained  for 
several  trial  runs.  Then  approximate  algebraic  expressions  of  the 
same  general  shape  were  assumed,  and  improved  by  trial  and  error 
until  the  estimation  was  judged  to  be  satisfactory.  In  this  manner, 
the  following  algebraic  equations  were  finally  obtained. 


Pn(t)  =  0,0025 


(3-41) 


P12(t)  =  0.007 


(3-42) 


P13(t)  =  -0.04  exp(-0.0016t) 


(3-43) 


P22(t)  =  0.0023 


(3-44) 


?23 ( t )  =  exp(-0.0016t) 


(3-45) 


P33(t)  =  exp(0.003t) 


(4-51) 


With  the  above  equations  replacing  the  covariance  differential 
equations  in  Filter  I,  the  process  was  simulated  with  the  previously 
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used  initial  state  and  noise  level.  For  ease  of  reference,  the  above 
filter  with  algebraic  expressions  for  the  elements  of  P_(t)  will  be 
referred  to  as  Filter  III.  Figure  3-6  shows  the  estimation  of  E-j  by 
Filter  III.  It  is  found  that  the  Filter  III  estimation  of  E-j  is 
better  than  the  estimates  provided  by  Filter  I  and  Filter  II.  Figure 
3-7  compares  the  estimation  of  c  obtained  using  the  three  filter 
algorithms  when  E-j(O)  contains  a  -10%  error.  A  comparison  of  Figure 
3-4  and  3-7  indicates  that  a  good  estimate  of  c  occurred  when  a  good 
estimate  of  E-j  was  obtained,  as  would  be  expected.  Table  3-3  compares 

A 

the  effect  of  E-j(O)  on  the  final  estimates  of  E,  c  and  T.  Even  with 
20%  errors  in  E-j(O),  the  estimation  of  E-j  by  the  simplified  algorithms, 
Filters  II  and  III  are  still  quite  encouraging.  The  poor  estimation 

A. 

of  extent  c  by  the  Filter  II  algorithm  when  E-j(O)  contains  a  +20% 
error  is  due  to  slower  convergent  rate  of  the  estimated  E-j ,  since 
when  E-j  is  higher  the  reaction  is  proceeding  slower,  and  hence  the 
estimated  c  is  lower  than  the  actual  extent. 

From  the  above  results,  the  simplified  nonlinear  filter  esti¬ 
mation,  seems  to  be  a  feasible  approach,  particularly  since  the  full 
filter,  Filter  I,  is  "suboptimal"  to  begin  with  due  to  approximations 
made  during  its  derivation. 

3.3.2  Effect  of  Update  Interval,  t  ,  on  Filter  II  Estimation 

In  the  above  simulation,  the  update  interval  for  the  covariance 
differential  equation,  t  =  10,  was  chosen  quite  arbitrarily.  In 

r 

actual  applications,  large  values  of  t  are  desirable  to  reduce 
computational  requi rements .  However,  if  t  is  too  large  the  estimation 
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Figure  3-6  Filter  III  Estimation  of  E-j  with  -  10%  Error  in  E,(0)(t  =  1  sec 
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Figure  3-7  Estimation  of  Extent  c  by  Three  Filters  with  -10%  Error  in  E-.  ( 0 ) ( t  =  1  sec 
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True  value  for  E-,  =  24000  cal/g-mole:  c ( 1 695 )  =  0.419,  T ( 1 695 )  =  666°K 
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may  not  converge.  Table  3-4  compares  the  effect  of  t  on  the  Filter 
II  estimates  at  the  switching  time  of  t  =  1695  sec.  Even  with  t  =  100, 

3  p 

the  Filter  II  algorithm  still  provides  a  good  parameter  estimate 


A 

although  the  estimated  extent,  c(1695),  deviates  farther  from  the 
actual  value  of  c (1 695 )  =  0.419  as  t  increases.  This  is  due  to  the 
fact  that  when  t  is  increased,  the  resulting  slower  convergence 
rate  of  estimated  E-j  produces  a  poor  estimate  of  extent  c.  A  further 
check  calculation  with  t  =  i.e.  P_(t)  =  P_(0)  for  all  t  >  0,  showed 
that  the  estimate  of  E-j  did  not  change  greatly  from  the  assumed  value 

A 

of  E-|(0)  and  the  estimated  extent  c  is  lower  than  the  actual  extent, 
as  would  be  expected. 


TABLE  3-4 


EFFECT  OF 

t  ON  PERFORMANCE  OF  FILTER  II 

P  .  .  

o 

o 

LQ 

CV1 

II 

o 

<UJ 

cal/g-mole, 

t  =  1  sec. ,  a  = 

1CC) 

E ( 1 695 ) (cal/g-mole) 

c  ( 1 695 ) * 

T (1695) (°K)* 

tp(sec 

24,043 

0.404 

667 

10 

23,909 

0.384 

668 

50 

23,988 

0.367 

668 

100 

^Actual  values  for  E-j  =  24000  cal/g-mole: 

c(l 695)  =  0.419,  T ( 1 695 )  =  666°K 
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3.3.3  Comparison  of  Computation  Times 

To  evaluate  how  much  computation  time  can  be  saved  by  using  a 
simplified  filter  algorithm,  object  decks  for  the  Fortran  program 
were  generated  for  all  three  filter  algorithms  and  estimates  of 
computation  times  were  obtained  for  an  IBM  360/67  computer  with  the 
MTS  operating  system.  Table  3-5  gives  an  appropriate  comparison, 
where  the  indicated  time  is  the  CPU  time  for  the  computer  program. 

It  is  found  that  the  CPU  time  can  be  reduced  by  58%  by  using  Filter 
II  with  t  =  10,  and  by  60%  using  Filter  III  in  comparison  with  the 
Filter  I  estimation.  Since  CPU  time  includes  both  program  I/O  and 
computation  times,  it  is  difficult  to  estimate  what  percentage  of 
CPU  time  is  concerned  with  the  computation  time. 

TABLE  3-5 

COMPARISON  OF  COMPUTATION  TIMES  FOR  THE  FILTER  ALGORITHMS* 

(E-j(O)  =  26400  cal/g-mole,  t$  =  1  sec.,  a  =  1°C) 

Filter  Computation  Time  (sec)  tp(sec) 


I 

49.3 

- 

II 

20.7 

10 

II 

18.1 

50 

II 

18,0 

100 

III 

19.6 

- 

This  value  includes 

1/0  operations. 

’ 

■ 
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3.3.4  Effect  of  Initial  State  Estimate 

In  some  physical  situations,  the  initial  state  is  not  known 
exactly.  For  example,  the  properties  of  the  initial  charge  to  a 
batch  reactor  may  change  in  an  unanticipated  manner,  particularly 
when  the  reactants  are  supplied  from  other  units  in  the  plant.  To 
see  how  an  inaccurate  initial  state  estimate  affects  the  filter  per¬ 
formance,  the  reaction  example  is  simulated  with  initial  estimates 
of  6(0)  =  0.01,  ^(0)  =  26400  and  two  values  of  f  (0) ,  490°K  and  510°K. 
In  this  simulation  the  previously  used  initial  covariance  matrix 
P_(0 )  was  used  for  both  Filter  I  and  II  algorithms.  Table  3-6  lists 
the  estimated  states  at  the  switching  time  while  Figure  3-8  shows  the 
response  of  the  estimated  activation  energy,  E-j(t).  The  response  is 
consistent  with  physical  intuition.  At  the  beginning  of  the  reaction 
c (0 )  =  0  and  T (0 )  =  500°K  but  the  estimated  values  are  c (0 )  =  0.01 
and  T(0)  =  510°K;  consequently,  the  estimated  E-j  increases  initially 
to  make  the  estimated  c(t)  and  T(t)  become  closer  to  the  actual 

A  A 

extent  and  temperature.  Finally,  when  c(t)  and  T(t)  approach  the 
actual  values  of  c(t)  and  T(t),  the  estimated  activation  energy, 

E(t),  converges  to  the  actual  value.  As  before,  better  parameter 
estimates  are  obtained  by  using  simplified  algorithms,  especially 
Filter  II  in  this  case. 
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TABLE  3-6 


(E1(0) 

Filter 

EFFECT  OF 

INCORRECT  INITIAL  STATE 

ESTIMATES 

,  0  =  1°C) 

T(1695)*(°K) 

ON  FILTER  PERFORMANCE 

=  26,400  cal/g-mole,  c(0)  =  0.01, 

f(0)(°K)  E1 (1695) (cal/g-mole) 

t  =  1  sec. 
s 

8(1695)* 

I 

510 

23,385 

0.453 

685 

I 

490 

23,331 

0.454 

686 

II 

510 

24,030 

0.408 

667 

II 

490 

24,011 

0.417 

667 

III 

510 

24,016 

0.444 

668 

III 

490 

23,923 

0.453 

666 

★ 

True  value  for  E-|  =  24,000  cal/g-mole 

c Cl 695 )  =  0.419  ,  T (1695)  =  666°K 

3.3.5  Effect  of  Noise  Level 

As  has  been  discussed  in  Section  3.3.1,  the  method  developed 
by  Chien  and  Aris  [5]  for  parameter  estimation  is  quite  sensitive  to 
noise  in  the  temperature  measurements.  Even  with  a  =  1°C,  they  did 
not  obtain  satisfactory  parameter  estimates.  To  see  how  large  a  noise 
level  the  nonlinear  filter  can  allow,  the  effect  of  noise  level  on 
Filters  I,  II  and  III  was  investigated.  From  simulation  results,  it 
was  found  that  a  large  noise  level  will  result  in  a  slowly  convergent 
estimation.  It  was  also  found  that  the  simplified  filter  algorithms 
can  handle  large  measurement  errors.  Table  3-7  shows  the  final 
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estimates  of  E,  c  and  T  by  Filters  I,  II  and  III. 

/V 

When  E-j(O)  has  -10%  error  and  ts  =  1 ,  the  Filter  I  estimation 
scheme  "blew  up"  because  of  numerical  instability  in  the  filter 

A 

equations.  When  the  estimated  activation  energy,  E-j  ,  is  too  low,  the 
fixed  integration  step  size  of  1  sec.  gave  large  computation  errors 
especially  in  the  evaluation  of  the  Jacobian.  As  a  further  check, 
calculations  were  made  by  using  a  data  sampling  time  of  t  =  0.5  sec. 

A 

for  the  Filter  I  estimation  when  E-j(O)  has  a  -10%  error  involved  and 
a  =  10°C.  As  indicated  in  Table  3-7,  for  t  =  0.5  the  numerical 
instability  was  eliminated  and  the  estimation  converged. 

From  a  practical  point  of  view,  a  =  10°C  represents  a  large 
temperature  measurement  error  and  so  higher  noise  levels  were  not 
investi gated. 

3.3.6  Effect  of  Data  Sampling  Interval 

In  all  the  previous  calculations,  a  data  sampling  time  of 

t  =  1  was  assumed  for  the  temperature  measurements.  In  real-time 
s 

applications,  large  sampling  intervals  are  desirable  to  reduce  the 
loading  on  the  computer.  However,  if  the  sampling  interval  is  too 
large,  the  filter  may  yield  poor  estimates.  To  test  the  "robustness" 
of  the  filter  algorithms,  a  data  sampling  time  of  tg  =  10  was 

A  A 

investigated.  The  standard  set  of  conditions  of  c(0)  =  0,  T(0)  =  500, 
a  =  1°C  and  the  same  P_(0 )  were  used  as  before.  Figures  3-9  and  3-10 
show  the  estimated  responses  with  -10%  and  +10%  errors  involved  in 
E  (0),  respectively.  When  compared  with  the  filter  estimation  for 
t  =  1  as  shown  in  Figures  3-4,  3-5  and  3-6,  the  convergence  rates 


. 


EFFECT  OF  NOISE  ON  FILTER  PERFORMANCE 
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of  E-|  estimation  are  essentially  the  same,  although  large  computation 
errors  are  involved  in  the  filter  calculations  due  to  the  large 
integration  step.  A  further  calculation  with  an  error  in  E-j(O)  of 
-10%  and  t  =  15  was  checked.  Both  the  Filter  I  and  II  estimates 
become  unstable,  while  Filter  III  still  worked  quite  well.  Since 
more  differential  equations  are  integrated  in  Filters  I  and  II  than 
in  Filter  III,  numerical  instabilities  appear  to  be  more  critical  in 
the  Filter  I  and  II  algorithms  when  the  integration  step  for  the 
filter  equations  is  increased  (i.e.  the  data  sampling  interval 
i ncreases ) . 

3.4  Concl usions 

(1)  The  optimal  control  policy  is  very  sensitive  to  parameter 
variation.  Hence  if  a  parameter  does  vary  and  is  undetected,  the 
optimal  control  policy  is  not  feasible. 

(2)  When  the  activation  energy  E-j  changed  and  was  not  iden¬ 
tified,  Mi  liman's  proportional  control  law  still  seems  to  be  practical, 
providing  the  extent  can  be  measured. 

(3)  In  using  Detchmendy  and  Sridhar's  nonlinear  filter, 
numerical  experimentation  is  required  to  find  a  satisfactory  P_(0 )  • 

In  finding  a  satisfactory  ?_{ 0)  for  the  single  reaction  example, 
the  diagonal  element  corresponding  to  the  parameter  is  much  larger 
than  the  elements  corresponding  to  the  states. 

(4)  Simplified  filter  algorithms  can  be  applied  to  estimate 
unknown  E-j  during  the  adiabatic  path.  It  was  found  the  simplified 
nonlinear  algorithms  are  less  sensitive  to  noise  level,  data  sampling 


>r  ( s ; 


II 
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time  and  error  in  the  initial  state  estimate. 

(5)  To  find  a  set  of  satisfactory  algebraic  equations  for  the 
Filter  III  algorithm,  one  approach  is  to  find  the  response  of  all 
elements  of  P_(t )  for  a  Filter  I  or  II  simulation,  then  exponential 
expressions  can  be  used  as  approximations  to  these  responses.  Since 
the  response  of  P_(t)  dependson  different  initial  estimates  of  states 
(including  the  initial  estimate  of  unknown  parameter),  the  responses 
of  P_(t)  are  required  for  several  initial  state  estimates  in  order  to 
obtain  satisfactory  algebraic  equations  for  Filter  III. 
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CHAPTER  IV 

STUDIES  OF  A  COMPLEX  REACTION 


4.1  Formulation 

Consider  the  complex  batch  reaction  originally  studied  by 
Rosenbrock  and  Storey  [10]  and  later  investigated  by  Millman  [12]: 


^4  ^4 

where  A-j  is  the  feed  material,  A^  an  intermediate,  A^  the  desired 
product  and  A^  the  wasteful  side  product. 

The  rate  constants,  k. ,  are  given  by 

k.  =  k!  tb  exp[-E.(l/T  -  1/658)]  (4-2) 

i  =  1 ,2 , ...  ,5 

where  k!  is  the  frequency  factor,  tb  the  total  batch  processing  time, 

E.  the  activation  energy  (divided  by  the  gas  constant)  and  T  is  the 
temperature . 

For  the  kinetic  parameters  given  in  Table  4-1,  the  optimal 
temperature  control  policy  to  give  a  maximum  concentration  of  A^  has 
been  studied  when  the  reaction  takes  place  in  a  plug-flow  tubular 
reactor  [10].  Millman  considered  batch  reactor  operation  and  com¬ 
pared  the  temperature  schedules  of  optimal  and  proportional  control 
policies  for  a  fixed  batch  processing  time  tb>  Following  Millman  [12], 
the  rate  equations  governing  the  state  of  the  reaction  are 


. 
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dx  -I 

rl  "  dt  " 

-(k]  +  k2  +  k3)  X-, 

(4-3) 

dx2 

r2  dt  " 

^4 0  ”X-!  -x2-x2 )  “  ^5xi 

(4-4) 

dx~ 

r3  dt  “ 

(k2  k^)  x-j  +  k^x2 

(4-5) 

wi  th 

0  - 

t  <  1 

where  time  t  is  dimensionless  due  to  normalization  with  the  total 
batch  processing  time  t^ ,  and  x-|  ,  x2  anc*  x3  are  concentrations  in 
moles  per  unit  volume  of  A-j  ,  Ag  and  A^  respectively.  The  concentra¬ 
tion  of  A2  can  be  determined  from  the  overall  mass  balance.  The 
state  variables  are  chosen  to  be  A-j  ,  A^  and  A^  since  it  is  assumed 
in  the  numerical  computations  of  Section  4.3  that  noisy  measurements 
of  these  variables  are  available. 

TABLE  4-1 

KINETIC  PARAMETERS  OF  THE  COMPLEX  REACTION 


Ei  (°K) 


.077 

8,000 

.070 

7,000 

.029 

7,500 

.248 

5,000 

.006 

7,500 

5 
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4.2  Control  Studies 

4.2.1  Suboptimal  Control  Policy 

It  has  been  discussed  in  Chapter  III  that  the  optimal  tempera¬ 
ture  control  policy  for  a  single  reversible  exothermic  reaction  tak¬ 
ing  place  in  a  batch  reactor  is  to  have  a  maximum  reaction  rate  at 
every  instant.  This  situation  will  not,  in  general,  result  in  a 
maximum  yield  of  the  desired  product  for  complex  reactions  [11]  and 
hence  is  not  an  optimal  policy.  Computational  techniques  such  as 
Pontryagin's  maximum  principle,  dynamic  programming,  hill  climbing 
methods,  etc.  can  be  applied  to  calculate  the  optimal  solution, 
providing  it  exists.  However,  since  the  required  calculations  are 
normally  quite  extensive  and  result  in  an  "open-loop"  control  policy, 
updating  the  optimal  control  on-line  to  compensate  for  process 
parameter  changes  is  generally  not  feasible.  Consequently,  only 
suboptimal  control  policies  which  can  be  adapted  on-line  will  be 
investigated. 

Two  suboptimal  control  policies  are  discussed  in  this  complex 
reaction  example;  one  will  be  referred  to  as  Mi  liman's  proportional 
control  law  and  the  other  as  the  instantaneous  performance  index 
control  law. 

With  all  the  kinetic  parameters  assumed  to  be  known  exactly 
as  in  Table  4-1,  Millman  [12]  investigated  a  suboptimal  control 
policy  consisting  of  the  following  proportional  control  law: 

T  =  Ki  +  K2  x2  (4_6) 


■  '  ■  '  V, 
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and  calculated  the  best  values  of  K-j  =  1120  and  K2  =  513  through  the 
use  of  a  gradient  method.  (This  control  law  was  also  considered  in 
Section  3.2.3  for  the  single  reaction.)  His  results  [12]  indicate 
that  this  suboptimal  control  resulted  in  yields  which  were  very  close 
to  the  optimal  yields  for  several  reaction  examples. 

A  second  type  of  suboptimal  control  policy  can  be  derived 
by  maximizing  a  performance  index  which  depends  only  on  the  current 
values  of  the  state  variables.  Several  instantaneous  performance 
index  control  schemes  have  been  discussed  by  Weber  and  Lapidus  [15]. 
Following  Seinfeld  [14]  a  suboptimal  control  policy  will  be  derived 
by  maximizing  the  instantaneous  rate  of  reaction.  The  resulting 
control  policy  will  be  referred  to  here  as  the  instantaneous  per¬ 
formance  index  ( I P I )  control  law. 

Since  is  the  desired  product  in  this  example,  the  suboptimal 
control  policy  based  on  maximizing  the  instantaneous  rate  of  reaction 
of  A^  will  be  the  solution  of 

ar0 

ar  -  0  (4'7) 

Thus  differentiating  Equation  (4-4)  with  respect  to  T,  the 
temperature  schedule  is  obtained  as 

T(t)  -  V<1n[|  (1.x^rx73/(E5-E4)  »  V658}  (4-8) 

A  comparison  of  the  suboptimal  and  optimal  control  policies 
performing  under  "ideal  conditions"  is  shown  in  Table  4-2.  In  the 
digital  computer  simulation  it  was  assumed  that  the  rate  constants 


' 


COMPARISON  OF  TEMPERATURE  CONTROL  POLICIES  FOR  THE  COMPLEX  REACTION 
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were  known  and  that  noise  free  measurements  of  x-j  ,  X2  and  x^  were 
available  at  a  data  sampling  time,  t  ,  of  0.005.  As  in  Millman's 
calculation  of  the  optimal  control  policy  [12],  the  temperature  was 
constrained  by,  850°K  <  T  <  1120°K.  The  results  in  Table  4-2 
indicate  that  the  suboptimal  policies  produce  final  concentrations 
of  which  are  very  close  to  the  optimal  value  of  0.448, 

4.2.2  Effect  of  Parameter  Sensitivity  on  the  Yield  of  A^ 

Again  it  is  assumed  that  the  initial  charge  characteristics 
change  unpredi ctably  from  batch  to  batch  and  can  be  represented  by 
random  variations  in  the  activation  energies,  and  Eg.  It  is  of 

Interest  to  determine  the  effect  of  undetected  changes  in  E^  and 
Eg  on  the  final  yield  of  Ag. 

In  Table  4-3  yields  are  compared  from  the  three  control 
strategies  when  the  actual  kinetic  parameters  have  changed  by  20% 
from  the  nominal  values  in  Table  4-1.  The  temperature  control 
policies  in  Table  4-3  are  based  on  the  nominal  parameter  values. 

The  "optimal  control  policy"  in  Table  4-3  was  approximated  by  linear 
interpolation  of  the  optimal  temperature  sequence  calculated  by 
Millman  [12]  and  shown  in  Table  4-2. 

It  is  found  that  the  optimal  temperature  schedule  based  on 
the  assumed  activation  energies  does  not  always  result  in  a  higher 
yield  than  suboptimal  control  policies,  when  the  activation 
energies  change  as  in  cases  1  and  2  of  Table  4-3.  One  disastrous 
situation  for  the  IPI  control  policy  is  case  3  where  a  yield  of  only 


0.174  was  obtained. 
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The  effect  of  variations  in  the  activation  energies  and 
on  the  control  schemes  has  just  been  discussed.  A  similar  situation 
occurs  when  inaccurate  parameter  values  are  used  in  the  calculation 
of  the  control  policy,  i.e.  suppose  that  the  process  parameters 
remain  at  the  nominal  values  of  Table  4-1,  but  the  assumed  model 
parameters  differ.  Table  4-4  shows  the  result  of  applying  the  IPI 
control  law  in  this  situation.  Only  in  case  4,  does  the  final  yield 
differ  significantly  from  the  optimal  yield  of  0.448  or  even  the 
suboptimal  yield  of  0.430  obtained  from  the  IPI  control  law  when  the 
activation  energies  are  known  exactly  (cf.  Table  4-2). 

This  type  of  comparison  was  not  made  for  the  proportional 
control  and  optimal  control  policies  on  account  of  the  extensive 
programming  effort  that  would  be  required  to  prepare  the  necessary 
computer  programs. 

4.3  Sequential  State  and  Parameter  Estimation 
4.3.1  Open  Loop  Parameter  Estimation 

Suboptimal  control  policies  with  noise  free  measurements  have 
been  investigated  in  the  previous  section.  It  was  found  that 
Millman's  proportional  control  law  is  not  very  sensitive  to  the 
variations  in  E4  and  E5-  However,  the  required  computations  to 
find  the  best  controller  constants,  K-j  and  K2,  for  a  given  set  of 
activation  energies  are  extensive.  The  IPI  control  policy  is  quite 
sensitive  when  the  assumed  values  of  E^  and  E&  used  in  the  control 
calculations  are  too  high  and  too  low,  respectively  (as  shown  in 
cases  3  and  4  of  Tables  4-3  and  4-4,  respectively). 
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Hence  on-line  parameter  estimation  may  be  desirable  to  improve  the 
final  yield  of  by  compensating  for  uncertain  or  changing  para¬ 
meters  . 

Assume  no  dynamic  error  is  involved  in  the  process  and  and 
are  the  only  undetermined  parameters.  Again  Detchmendy  and 
Sridhar  nonlinear  filter  technique  [22]  is  applied  to  perform  the 
on-line  parameter  (and  state)  estimation. 

It  is  assumed  that  the  measured  values  of  the  concentrations 
of  A-j  ,  A^  and  A^  contain  random  noise  as  described  by  the  error  model: 


y.  =  x -  [1  +  e_*  (0 ,a) ]  ,  i  =  l  ,2,3 


(4-9) 


where  y^  is  the  measured  concentration  of  x^.  and  the  noise  error  e,  is 
assumed  to  be  Gaussian  distributed  with  zero  mean  and  a  standard 
deviation,  a. 


Define  state  variables  x-j  (t) ,  x2(t)  and  x3(t)  as  the  corres¬ 


ponding  estimated  values  of  x-j  ,  x^  and  x^  at  time  t.  To  reduce  the 
"stiffness  effect"  in  the  filter  equations  (as  in  the  simple 
reaction  example),  state  variables  corresponding  to  E4  and  E5  are 


defined  as  x^(t)  and  x^(t)  with 


E4(t)  -  E4,ref 

E4,ref 


(4-10) 


i5(t) 


5  ,ref 


(4-11) 


where  E.(t)  and  E5(t)  are  the  estimated  values  of  E4>  E&  at  any 


. 

. 
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instant  t  and  re^s  ^  are  fixed  reference  values. 

Then  the  filter  equations  (2-7  to  2-9)  can  be  developed  with 
the  state  equations  and  Jacobian  as 


A 


— 

—  ( k  -|  +  ^2  kg)  X-J 

(4-12) 

f2 

= 

/V  AAA  A  A 

*<4(1  xi  x2  x3^  “  ^5  x2 

(4-13) 

f3 

= 

A  A*  A  A 

(l<2  +  )  x-j  +  kg  Xg 

(4-14) 

f4 

= 

0 

(4-15) 

f5 

= 

0 

(4-16) 

fll 

- 

—  ( k  1  +  k£  -*■  kg) 

(4-17) 

f12 

= 

0 

(4-18) 

fl  3 

- 

0 

(4-19) 

f14 

= 

0 

(4-20) 

f15 

= 

0 

(4-21) 

f21 

- 

A 

-  k4 

(4-22) 

f22 

= 

;a  a 

"  k4  -  k5 

(4-23) 

f23 

= 

A 

-  k4 

(4-24) 

f24 

- 

(l-^-Xg-Xg)  k4  [-  E4,refU/T  '  1/658)] 

(4-25) 

f25 

= 

~*2  h  C-  E5,ref  <1/T  '  1/658)] 

(4-26) 
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f31 

- 

k2  +  £ 

4 

(4-27) 

f32 

= 

A 

k5 

(4-28) 

f33 

- 

0 

(4-29) 

f34 

= 

*3" 

<  -X 

<  X 

E4,  ref 

(1/T 

-  1/658)] 

(4-30) 

f35 

= 

A  A 

x2  k5 

^  ^5, ref 

(1/T 

-  1/658)] 

(4-31) 

f41 

= 

f42  = 

f43  =  f44 

=  f45 

=  0 

(4-32) 

f51 

= 

f52  = 

f53  =  f54 

=  f55 

=  0 

(4-33) 

with 

the  re‘ 

lati 

ons 

A 

k4 

= 

tb) 

exp  [ 

^E4,ref 

A 

+  x4 

W**1'1  - 

1/658)] 

(4-34) 

A 

k5 

=  (k5 

tb) 

exp  [ 

'^5 ,ref 

A 

+  x5 

E5.refW/T  - 

1/658)] 

(4-35) 

Suppose  the  reactor  state  is  known  initially  and  the  reference 
values  of  ^  and  Eg  are  chosen  to  be  the  initially  assumed 
values  of  E^  and  Eg.  Then  the  initial  conditions  of  the  filter 
equations  are, 

^(0)  =  1 

x2(0)  -  x3(0)  =  x4(0)  =  x5(0)  -  0 

irrespective  of  the  assumed  values  of  E^  and  Eg. 

For  this  investigation  of  the  nonlinear  filters  under  open 
loop  conditions,  it  will  be  assumed  that  the  reactor  temperature 


• 

. 
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varies  linearly  with  time: 

(obtained  from  a  linear  approximation  to  the  optimal  temperature 
schedule  of  Table  4-2) 


T  =  1120  -  250  t 


(4-36) 


As  in  Section  3.3.1  numerical  experimentation  was  required  to 
find  a  satisfactory  initial  covariance  matrix  P_(0) .  From  the 
experience  of  finding  P_(0)  in  Section  3.3.1  for  the  single  reactor, 
a  diagonal  matrix  of  £(0)  was  first  tested  with  the  elements  cor¬ 
responding  to  parameters  larger  than  the  elements  corresponding  to 
the  state  variables.  A  slow  convergent  rate  of  filter  estimation 
resulted  using  a  diagonal  matrices  for  £(0).  To  improve  the  rate  of 
convergence,  the  off-diagonal  elements  of  £(0)  were  also  tested. 
After  20  trial  runs,  the  following  initial  covariance  matrix  £(0) 
was  found  to  be  satisfactory. 


P(0) 


0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

0.5 

250 

oi 

0.5 

0.5 

0.5 

0.5 

250 

As  discussed  in  Section  2.2,  the  weighting  matrices  in  the 
performance  index  must  also  be  specified.  Since  dynamic  errors  are 
not  considered  in  this  complex  batch  reaction,  only  the  weighting 
matrix,  Q_(t),  must  be  specified.  For  this  example,  the  following 


matrix 
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• 
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act) 


1  0  0 
0  1  0 
0  0  1 


with  0  <  t  s  1 


was  arbitrarily  chosen  to  place  identical  weight  on  each  measured 
concentration. 

With  all  the  initial  conditions  available,  the  filter  equations 
were  integrated  by  a  Runge-Kutta-Gi 1 1  4-th  order  routine  [43]  with  the 
integration  interval  equal  to  the  data  sampling  interval  of  t  =  0.005. 
The  simulated  process  data  were  generated  by  integrating  Equations 
(4-3)  to  (4-5)  with  the  same  integration  routine  and  the  measured 
concentrations  calculated  from  Equation  (4-9)  using  unbiased  Gaussian 
distributed  random  numbers  with  a  standard  deviation  of  a  =  0.1  and 
zero  mean.  For  the  piecewise  recursive  filter  algorithm,  t  =  0.01 
was  used,  i.e.  the  updating  interval  of  the  covariance  differential 
matrix  was  twice  the  sampling  interval.  As  in  Chapter  III,  the 
original  filter  algorithm  will  be  referred  to  as  Filter  I  and  the 
piecewise  recursive  filter  algorithm  as  Filter  II.  The  filter  esti¬ 
mates  are  shown  in  Figures  4-1  and  4-2  for  20%  errors  in  the  initial 
guesses  of  E^  and  Eg.  The  solid  and  dashed  lines  correspond  to 
cases  1  and  2  respectively  of  Table  4-4. 

The  results  in  Figure  4-1  and  4-2  indicate  that  the  estimates 
of  E^  converge  quite  rapidly  for  both  filters  while  the  estimates  of 
Eg  converge  slowly  for  both  filters  in  case  2  where  the  initial  guess 
of  Eg  is  20%  low.  For  the  initial  guesses  of  E4  and  Eg  corresponding 
to  cases  3  and  4  of  Table  4-4,  the  results  are  quite  similar  to 
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Figure  4-2  Open  Loop  Estimation  of  E 
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Figures  4-1  and  4-2.  The  estimation  of  converges  very  quickly, 
while  the  Eg  estimation  converges  very  slowly  especially  when  -20% 
error  involved  in  initial  guesses  of  Eg.  The  reason  for  a  slowly 
convergent  estimate  of  Eg  is  due  to  the  reaction  mechanism.  When  the 
reaction  first  begins,  the  formation  of  from  the  reaction  path  of 
k2  and  k3  as  shown  in  Equation  (4-1)  dominates,  hence  the  estimated  ' 
concentration  of  does  not  differ  greatly  from  the  actual  concentra¬ 
tions  of  A4  in  the  beginning,  thus  the  driving  term,  [y_  -  h_(x,t)],  in 
the  filter  equations  is  very  small. 

In  Filter  I  estimation,  5  state  variables  are  estimated 
(including  the  two  unknown  parameters).  Therefore,  20  differential 
equations  are  integrated  every  time  a  sample  is  taken.  If  algebraic 
equations  can  be  used  to  replace  the  15  covariance  differential 
equations  of  Filter  I,  as  suggested  by  Seinfeld  and  Chen  [14],  a 
significant  saving  in  computation  time  is  predictable.  In  order  to 
obtain  satisfactory  algebraic  equations,  the  time  history  of  P_(t) 
from  a  Filter  I  or  Filter  II  simulation  is  needed,  then  exponential 
approximations  are  assumed  to  approximate  the  response  of  each 
element  of  P(t).  If  these  approximate  expressions  provide  poor  para¬ 
meter  estimates,  numerical  experimentation  is  again  required  to 
modify  the  algebraic  equations.  About  10  trials  were  required  to 
determine  the  following  algebraic  equations  which  were  chosen  to 
replace  the  covariance  differential  equations. 

Pn(t)  =  P12(t)  =  P13(t)  =  P14(t)  =  P15(t)  =  0  (4-36) 

P22(t)  =  19.8  exp(-3.3t) 


(4-37) 


■ 


. 


71 

P23(t)  =  "3*4  exPH  .5t) 

(4-38) 

P24(t)  =  117  exp(-8.8t) 

(4-39) 

P25U)  =  -17.8  exp(-l  .21 1) 

(4-40) 

P33 ( t )  =  9*8  exp(-2.1t) 

(4-41) 

P34U)  =  8.2  exp(-l  .44t) 

(4-42) 

P35(t)  =  47.6  exp ( - 1 . 85t ) 

(4-43) 

P44(t)  =  270  exp(-3.0t) 

(4-44) 

P45 ( t )  =  43.8  exp(-3.7t) 

(4-45) 

P55(t)  =  533  exp ( -4 . 7t ) 

(4-46) 

The  filter  obtained  by  replacing  the  elements  of  P_(t)  with 
algebraic  equations  will  be  referred  to  as  the  Filter  III  algorithms 
Figure  4-3  shows  the  estimation  provided  by  Filter  III  with  initial 
guesses  of  and  corresponding  to  cases  1  and  2  of  Table  4-4.  By 
comparison  with  the  Figures  4-1  and  4-2,  the  Filter  III  algorithm 
seems  to  be  a  feasible  approach  to  perform  unknown  parameter  esti¬ 
mation.  The  same  situation  in  the  fast  convergence  of  E^  is  quite 
predictable,  since  the  algebraic  equations  are  chosen  to  approximate 
the  time  history  of  the  covariance  matrix  elements  of  Filters  I  and  II. 

4.3.2  Comparison  of  Computation  Times  for  Filters  I,  II  and  III 

In  order  to  get  a  rough  idea  about  how  much  computation  time 
can  be  saved  by  using  the  simplified  nonlinear  filter  algorithms, 
object  decks  for  Filter  I,  II  and  III  were  generated  to  perform  the 
computation  on  an  IBM  360/67  computer  with  the  MTS  operating  system. 
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Table  4-5  gives  the  computer  CPU  time  of  Filter  I,  II  and  III  estima¬ 
tion  with  initial  estimates  of  and  corresponding  to  case  1  of 
Table  4-4  and  tg  =  0.005.  It  is  found  that  a  25%  reduction  in  computa¬ 
tion  time  can  be  obtained  by  using  Filter  II  with  t  =  2  t  ,  while  the 

p  s 

Filter  III  algorithm  can  save  40%  of  the  computation  time  in  compari¬ 
son  with  the  Filter  I  estimation. 


TABLE  4-5 

COMPARISON  OF  COMPUTATION  TIMES 


FOR  THE 

FILTER  ALGORITHMS  (t  =  0.005) 

Filter 

* 

Computation  Time  (sec) 

*p 

I 

21 .8 

II 

16.3 

2  t 

s 

II 

13.8 

4ts 

III 

12.2 

_ 

includes  I/O  operations 


4.3.3  Effect  of  Sampling  Interval  on  Parameter  Estimation 

Since  most  of  the  real-time  computers  used  in  the  process 
industries  are  small  or  medium  size  computers,  the  computational 
requirements  made  of  the  computer  is  another  important  aspect  of 
on-line  parameter  estimation.  If  a  large  sampling  time  can  be  used, 
the  load  on  the  computer  for  data  acquisition  and  the  filter  calcu¬ 
lations  will  be  reduced.  Hence  the  effect  of  sampling  interval  on 
parameter  estimation  will  be  investigated  for  all  three  filter 
algorithms.  For  t  =  0.01,  the  resulting  parameter  estimates  showed 
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essentially  the  same  convergence  as  indicated  in  Figures  4-1  to  4-3 
for  ts  =  0.005.  As  a  further  check,  t  =  0.02  was  investigated  for 
all  three  filters  with  t  =  2  t$  being  used  for  the  Filter  II  estima¬ 
tion.  Figure  4-4  shows  typical  results  for  the  case  1  initial  con¬ 
ditions  of  Table  4-4.  All  filters  provide  a  rapidly  converging 
estimate  of  E^,  while  the  estimation  of  is  poor  for  Filters  I  and 

III,  but  Filter  II  still  gives  a  reasonable  estimation.  When  t  was 

s 

increased  to  a  value  of  0.04,  all  three  filters  gave  poor  estimates 
of  E^  and  E^ . 

4.3.4  Effect  of  t  on  Filter  II  Estimation 

To  investigate  the  effect  of  t  on  Filter  II  estimation, 

r 

values  of  t  =  0.005  and  t  =  0.02  (i.e.  t  =  4  t  )  were  used  in  the 

5  p  P  S 

A  A 

simulation.  For  the  initial  parameter  estimates  of  E^(0)  and  Eg (0 ) 

corresponding  to  cases  1  and  2  of  Table  4-4,  Figure  4-5  shows  the 

Filter  II  estimation  with  t  =  4  t  .  It  is  found  that  the  parameter 

P  * 

estimates  are  poor  even  for  E^.  By  comparison,  in  the  last  section, 

when  t  =  0.02  and  t  =  0.04  (i.e.  t  =  2  t  )  were  used  in  the  Filter 
s  p  p  s 

II  estimation,  good  estimates  of  E^  and  E5  were  obtained.  Hence  in 
the  actual  application  of  the  Filter  II  algorithm,  the  relation  of 
t  =  2t$  appears  desirable  in  this  example  in  order  to  obtain  satis¬ 
factory  estimates.  Despite  this  restriction  on  t  ,  a  significant 
reduction  in  computer  time  is  achieved  by  Filter  II  as  shown  in 
Table  4-5. 

4.3.5  Parameter  Estimation  When  Some  States  are  not  Measured 

In  the  above  calculations,  the  state  knowledge  required  in 
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Figure  4-4  Open  Loop  Estimation  of  E4  and  Eg  by  the  Three  Filters  (case  1,  t  =  0.02 


case 
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Figure  4-5  Effect  of  t  on  Filter  II  Estimation  (t  =  0.005, 
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control  policy  is  assumed  to  be  measurable.  In  many  processes,  it  is 
either  not  possible  or  practical  to  measure  all  of  the  state  variables. 
In  this  situation,  the  nonlinear  filter  can  be  used  to  provide  esti¬ 
mates  of  all  the  states  from  measurement  data  of  the  states  which  are 
accessible.  To  evaluate  the  applicability  of  the  nonlinear  filter 
technique  when  some  states  are  inaccessible,  the  previous  complex 
reaction  scheme  is  examined.  Suppose  only  A-j  and  can  be  measured 
and  the  measured  states  follow  Equation  (4-4).  Let  and  E^  be  the 
undetermined  parameters.  Then  define  the  state  variables  and  develop 
the  filter  equations  (2-7  to  2-9)  as  in  Section  4.3.1  with  the 
arbitrarily  chosen  weighting  matrix 

0  S  t  S  1 

and  the  temperature  schedule  in  Equation  (4-36).  Then  the  process  is 
simulated  on  the  digital  computer  as  in  the  previous  computation 
scheme.  It  was  found  that  a  modification  of  the  previously  used  £(0) 
by  setting  P44(0)  =  50  and  P55(0)  =  75  resulted  in  essentially  the 
same  estimation  of  E-,  and  E5  as  in  case  1  of  Figure  4-1  by  Filter  I 
estimation.  But  when  the  initial  estimates  of  state  are  the  same  as 
case  2  of  Figure  4-1,  the  modified  P_(0)  gave  a  slowly  convergent 
estimation  of  E5  although  the  E4  estimation  converged  rapidly.  Further 
experimentation  by  adjusting  P44  and  P55  in  £(0)  could  probably  improve 
the  estimation. 

Hence  when  certain  states  are  required  to  implement  the  control 


Q(t)  = 
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policy  but  are  not  measured,  the  nonlinear  filter  technique  is  a 
feasible  way  to  supply  the  required  state  estimates. 

4.3.6  Closed  Loop  Parameter  Estimation 

In  Chapter  III,  the  parameter  estimation  and  optimal  cooling 
paths  could  be  separated  since  the  parameter  estimation  was  essentially 
completed  during  adiabatic  operation  (i.e.  before  the  switching  point 
for  optimal  cooling  was  reached).  In  this  complex  reaction  example, 
parameter  estimation  could  in  principle  be  performed  during  isothermal 
operation  or  some  other  specified  temperature  schedule  and  then  a 
switch  made  to  a  different  control  strategy  after  the  parameter  estima¬ 
tion  converges.  Unfortunately,  the  previous  open  loop  calculations 
indicated  that  the  estimation  of  inherently  involve  slow  convergence. 
An  alternative  control  strategy  proposed  by  Seinfeld  [14,38]  is  to  use 
feedback  control  with  a  nonlinear  filter  added  to  the  control  loop. 

In  this  case  the  parameter  estimation  and  control  action  are  performed 
simultaneously.  The  block  diagram  for  the  proposed  control  scheme  is 
shown  in  Figure  4-6.  Where  e(t)  is  the  observation  noise  and  the 
controller  can  be  any  suboptimal  controller.  Note  that  the  suboptimal 
control  is  based  on  the  estimate  of  the  state  supplied  by  the  filter 
rather  than  the  noisy  measurements. 

This  control  strategy  is  simulated  in  the  complex  reaction 
example  with  the  IPI  control  policy  chosen  as  the  suboptimal  control 
policy.  Assume  no  dynamic  errors  are  involved  in  the  process  and 
the  initial  state  of  the  reactor  is  known  exactly.  Then  with  the  same 
initial  covariance  matrix  P_(0)  and  weighting  matrix  (^(t)  that  were  used 
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Figure  4-6  Block  Diagram 
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in  the  open  loop  estimation,  the  closed  loop  system  is  simulated. 

Filter  I  estimations  are  shown  in  Figures  4-7  and  4-8,  Filter  II 
estimations  are  shown  in  Figures  4-9  and  4-10,  and  Figures  4-11  and 
4-12  show  the  Filter  III  estimations  with  initial  estimates  of 

A  /\ 

and  (i.e.  E^(0)  and  Eg (0 ) )  corresponding  to  the  4  cases  in  Table 

/V 

4-4.  As  in  the  previous  open  loop  estimation,  convergence  of  E^ 

A 

is  slower  than  that  of  E^.  In  comparison  with  the  previous  open  loop 
estimation,  Filter  I  and  III  provide  better  parameter  estimates 
during  close  loop  operation.  When  the  initial  estimate  of  E^  is 
low,  the  Filter  II  estimation  of  E^  is  poor.  The  Filter  II  estima¬ 
tion  could  be  improved  by  selecting  P_(0)  so  as  to  give  satisfactory 
parameter  estimates  by  Filter  II  during  closed  loop  operation.  As 
noted  earlier,  all  the  Filter  II  results  used  the  P_(0)  selected  for 
open  loop  estimation  by  Filter  I  (see  Section  4.3.1). 

Table  4-6  gives  the  final  yield  of  A^  when  the  nonlinear 
filter  is  added  to  the  control  loop.  In  comparison  with  the  results 
of  Table  4-4  where  a  nonlinear  filter  is  not  used,  the  addition  of 
the  nonlinear  filter  in  the  control  loop  only  small  changes  for  cases 
1  -  3.  But  when  the  initial  estimates  of  E4  and  E5  corresponded 
to  case  4  of  Table  4-4,  a  large  improvement  in  the  final  yield  of  A3 
was  obtained  when  the  nonlinear  filter  is  added  to  the  control  loop. 
Hence  when  the  parameters  are  not  known  exactly,  closed  loop  control 
with  a  nonlinear  filter  included  seems  to  be  practical. 

Another  problem  as  shown  in  Table  4-3  occurs  when  the  values 
of  E,  and  E,  change  (e.g.  due  to  catalyst  decay)  as  in  case  3.  Here 

H-  3 

the  IPI  control  policy  gives  a  very  poor  final  yield  of  A^  if  parameter 
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Figure  4-9  Closed  Loop  Estimation  of  E4  and  E5  by  Filter  II  for  Cases  1  and  2 
(t.  =  0.005,  a  =  0.1 ,  t  =  2t  ) 


ca  se  4 
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Figure  4-10  Closed  Loop  Estimation  of  E4  and  E5  by  Filter  II  for  Cases  3  and  4 
(tc  =  0.005,  a  =  0.1,  t  =  2t  ) 


case 
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Figure  4-11  Closed  Loop  Estimation  of  E4  and  Eg  by  Filter  III  for  Cases  1  and  2 
(t  =  0.005,  o  =  0.1) 


case 
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Figure  4-12  Closed  Loop  Estimation  of  E4  and  E-  by  Filter  III  for  Cases  3  and  4 
(t  =  0.005,  o  =  0.1) 
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TABLE  4-6 


YIELDS 

of  a3  OBTAINED  USING  CLOSED 

LOOP  PARAMETER 

ESTIMATION 

AND  THE  IPI 

CONTROL  POLICY 

(t  =  0.005, 

a  =  0.1) 

Percent  Error 

Percent  Error 

Final  Yield 

Estimated  Final 

Fi  1  ter 

in  E4(0) 

in  E5(0) 

of  A3 

Yield  of  A^ 

I 

+20% 

+20% 

0.437 

0.435 

I 

-20% 

-20% 

0,422 

0.426 

I 

-20% 

+20% 

0.431 

0.431 

I 

+15% 

-15% 

0.428 

0.431 

* 

II 

+20% 

+20% 

0.432 

0.431 

* 

II 

-20% 

-20% 

0.420 

0.421 

* 

II 

-20% 

+20% 

0.432 

0.432 

* 

II 

+15% 

-15% 

0.425 

0.427 

III 

+20% 

+20% 

0,432 

0.430 

III 

-20% 

-20% 

0.421 

0.426 

III 

-20% 

+20% 

0.434 

0.433 

III 

+15% 

-15% 

0.426 

0.429 

• 
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estimation  is  not  performed.  Table  4-7  compares  the  final  yields  of 
Ag  and  estimates  of  and  Eg  when  a  filter  is  added  and  when  no 
filter  is  used.  The  corresponding  concentration  histories  of  are 
shown  in  Figure  4-13.  The  estimated  concentration  of  A3  is  not  shown 
in  Figure  4-13,  since  the  estimated  values  are  quite  close  to  the 
actual  concentrations.  Figure  4-14  shows  the  Filter  I  estimation  of 
E^  and  Eg.  Thus  again  it  has  been  demonstrated  that  a  nonlinear 
filter  can  be  used  to  "track"  sensitive  parameters  during  batch 
reactor  operation  and  give  improved  control. 


TABLE  4-7 


Filter 

EFFECT  OF  VARIATIONS  IN  E4  AND  Eg  ON  FINAL  YIELD 

OF  A3  (IPI  CONTROL  POLICY,  t  =  0.005,  a  =  0.1) 

Final  Parameter 

Final  Yield  Estimated  Final  Percent  Error  in 
of  A3  Yield  of  A3  E4 

* k 

Estimates 

Percent  Error 
in  E5 

None 

0.173 

- 

-20% 

+20% 

I 

0.322 

0.325 

+2.4% 

+2.1% 

** 

II 

0.324 

0.319 

-1  .0% 

+9.6% 

III 

0.321 

0.322 

+1  .6% 

+5.0% 

★ 

Initial  parameter 

/V 

values  (i.e.  E^(0) 

estimates  were 

=  5000,  E5(0) 

chosen  to  be  the  nominal 

=  7500) 

parameter 

** 


. 

3 
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Figure  4-13  Yield  of  When  IPI  Control  Policy  is  Used  Both  With  Parameter 
Estimation  and  Without  Parameter  Estimation 
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Figure  4-14  IPI  Control  Scheme  with  Filter  I  Added  in  the  Feedback  Loop  to  Estimate 
E.  and  E-  (t„  =  0.005,  o  =  0.1) 
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4.4  Concl usions 

The  following  conclusions  can  be  made  for  the  complex  reaction 
example. 

(1)  The  parameter  sensitivity  analysis  showed  that  the  optimal 
and  proportional  control  policies  are  not  sensitive  to  parameter 
variations  while  the  IPI  control  policy  is  sensitive  to  only  one  com¬ 
bination  of  and  Eg  values  (case  3) .  It  also  showed  that  when  the 
parameters  varied  and  were  undetected,  the  optimal  control  policy  did 
not  always  result  in  a  higher  yield  of  desired  product  of  Ag  in  com¬ 
parison  with  the  suboptimal  control  policies. 

(2)  In  open  loop  parameter  estimation  by  Detchmendy  and 
Sridhar  nonlinear  filter,  numerical  experimentation  is  required 
to  find  a  satisfactory  P_(0).  It  was  found  that  the  diagonal  matrix 
JP(0)  resulted  in  slow  convergence  but  the  use  of  nonzero  off-diagonal 
elements  improved  convergence. 

(3)  For  some  initial  parameter  estimates,  the  estimation  of 

E  by  Filter  II  is  poor.  This  is  because  the  initial  matrix  f_(0)  used 

J 

in  Filter  II  was  based  on  numerical  experimentation  with  Filter  I .  In 
actual  applications,  numerical  experimentation  to  find  £_(0)  should  use 
Filter  II  directly. 

(4)  In  finding  a  satisfactory  algebraic  approximation  for  £.(0.) 
in  the  Filter  III  algorithm,  the  time  variation  of  P_(t)  elements  of 
Filter  I  or  II  are  needed  for  different  initial  estimates  of  state. 
Since  15  algebraic  equations  are  needed  in  the  example,  a  systematic 
way  of  finding  the  satisfactory  exponential  approximation  is  to 
determine  the  most  sensitive  elements  first,  then  numerical 


■ 

" 

. 


92 


experimentation  can  be  performed  on  the  sensitive  elements  to  find 
satisfactory  expressions. 

(5)  Closed  loop  scheme  of  performing  parameter  estimation 
and  control  action  simultaneously  seems  to  be  a  feasible  method  for 
"tracking"  changing  parameters. 
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CHAPTER  V 
CONCLUSIONS 

Since  most  physical  processes  are  subject  to  some  degree  of 
uncertainty  in  process  parameters,  a  parameter  sensitivity  analysis 
will  give  an  accuracy  requirement  for  the  control  policy.  The 
parameter  sensitivity  of  optimal  and  suboptimal  control  policies  of 
two  batch  reactions  have  been  presented  in  Chapters  III  and  IV 
respectively.  It  was  found  that  Mi  1 1  man  1 s  proportional  control 
policy  is  less  sensitive  than  the  optimal  control  policy  to  changes 
in  kinetic  parameters.  For  the  complex  reaction  example,  control 
policy  based  on  an  instantaneous  performance  index  is  quite  sensitive 
for  certain  conditions.  It  was  found  that  the  suboptimal  control 
policies  result  in  final  yields  which  are  quite  close  to  the  optimal 
control  policy  when  no  noise  is  involved  in  the  measu regents  and  the 
kinetic  parameters  are  known  a  priori.  Hence  suboptimal  control 
policies,  as  investigated  in  this  thesis,  seem  to  be  a  practical 
approach . 

For  the  batch  reaction  examples  investigated  in  this  thesis, 
the  proposed  nonlinear  filter  algorithms  provide  a  feasible  approach 
for  tracking  unknown  parameters.  Even  for  large  sampling  time  inter¬ 
vals,  high  measurement  noise  levels  and  inaccurate  initial  states, 
the  filter  algorithms  resulted  in  encouraging  parameter  estimations. 
When  the  computation  time  and  computer  storage  are  critical  factors, 
the  simplified  filter  algorithms  provide  a  reasonable  approach. 

A  filter  included  in  the  feedback  control  loop  can  track 
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uncertain  or  changing  parameters  and  improve  the  final  yields  of  the 
desired  products. 

Since  the  suboptimal  control  policies  require  knowledge  of  the 
current  state  of  the  process,  the  control  policies  cannot  be  imple¬ 
mented  when  some  of  the  required  state  variables  are  inaccessible. 
However,  a  nonlinear  filter  in  the  control  loop  can  also  provide  the 
required  state  estimates  from  the  measurement  data  and  thus  overcome 
this  shortcoming  of  many  suboptimal  control  policies. 
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NOTATION 

Symbol 

Description 

Ai 

i-th  chemical  species 

c 

extent  of  reaction 

c. ,  c. 

concentration  of  A.. 

cp 

Ei 

A 

Ei 

total  heat  capacity  per  unit  volume 

activation  energy 

estimated  E.. 

Ei  ,ref 
f(x.t) 

fixed  reference  value  of  E. 

vector  function 

fi 

i-th  component  of  f_(x_,t) 

fu 

9f.j  (x,t)/3Xj 

£(z,£st) 

vector  function 

h.(x,t) 

vector  function 

(AH) 

heat  of  reaction  (negative  for  exothermic  reaction) 

J 

(-AH)/Cp 

k  ' 

Ki 

frequency  factor 

K(x,t) 

matrix  function 

K1  ,K2 

I 

proportional  controller  constant 

dimension  of  unknown  input  vector 

m 

dimension  of  measurement  vector 

n 

dimension  of  state  vector  (including  unknown 
parameter) 

£ 

parameter  vector 

p(t) 

covariance  matrix  of  state  vector  x(t) 

P.  . 
ij 

element  of  P_( t ) 

q' 

rate  of  heat  removal 
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q 

fi(t) 

r(c,T),r. 

rm(c’T) 

R 

t 

lb 

tf 

S 

‘s 

T 

a 

T 

T 

m 

v(t) 

i(x.t) 

w(t) 

l(t) 

x(t) 

xi(t) 

X(t) 

x^t) 

*(t) 

z.(t) 


optimal  value  of  q' 
q,/cP 

optimal  value  of  q 
weighting  matrix 
reaction  rates 
optimal  value  of  r(c,T) 
gas  constant 
time 

total  batch  processing  time 
fixed  final  time 
update  interval  of  P_[t) 
data  sampling  interval 
temperature 
estimated  temperature 
optimal  value  of  T 
measurement  noise  vector 
matrix  function 
unknown  input  vector 
weighting  matrix 

argumented  state  vector  (including  parameters) 
i -th  component  of  x(t) 
estimate  of  x^(t) 

A 

i-th  component  of  x_(t) 
measurement  vector 
state  vector 
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0. 

l 

Yi 

P 

a 

e(0,a),e.j (0, a) 


stoichiometric  coefficient 
forward  reaction  order 
backward  reaction  order 

-  Vi 

standard  deviation 

Gaussian  distributed  random  noise  with  zero  mean 
and  standard  deviation  a 
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APPENDIX  A 

DERIVATION  OF  DETCHMENDY  AND  SRIDHAR  NONLINEAR  FILTER  [22,31] 

In  this  appendix,  Detchmendy  and  Sridhar  nonlinear  filter  [22] 
based  on  the  derivation  presented  by  Sage  [31]  is  derived.  First,  the 
Pontryagin  maximum  principle  corresponding  to  the  Bolza  problem  in 
the  calculus  of  variations  is  discussed.  Then  the  optimal  estimation 
problem  using  a  weighted  least  squares  criterion  is  formulated. 

Finally  the  invariant  imbedding  technique  is  employed  to  solve  the 
resulting  two  point  boundary  value  problem.  After  certain  approxi¬ 
mations  are  made,  the  final  filter  equations  are  resulted. 


The  Bolza  Problem  In  The  Calculus  Of  Variations 

Consider  a  given  nonlinear  system  operating  over  the  fixed 
time  interval  t  <  t  <  tf  of  the  form 

x  =  f(x,u,t)  (A_1) 


where  x(t) ,  the  n-vector  of  state  variables,  is  determined  by  u.(t), 

the  m-vector  of  control  variables. 

It  is  desired  to  determine  the  control  vector  u_( t )  such  that 

the  following  cost  function  is  minimized 


I 


e[x(t),t] 


d>[x(t),  u(t),  t]  dt 


(A-2) 


Next,  the  Lagrange  multiplier  _L_is  used  to  adjoin  the  system 
differential  equality  constraint  (i.e.  Equation  (A-l))  to  the  cost 


function,  which  gives 


t;  )u.y.d  b9nffin93-3b  z\  ^asfds'nsv  ^622  to  W:>$v*rr-IH&  <t)x  stdflw 
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I  =  e[x,(t),t] 


t=t. 


t=to  ~° 


U[x(t)  ,u(t)  ,t]+A '  [f  (x,t)-xj}dt  (A-3) 


where  X1  is  the  transpose  of  A_  . 

Define  a  scalar  Hamiltonian  function  as 

H[x.  =  <j>[x(t),u_(t),t]  +  l'[f(x,u.,t)-x] 

Then  the  cost  function  becomes 


(A-4) 


I  =  e[x,t] 


t=t. 


t=t  J  t 

O  0 


{H[x  ,u , A ,t]-A 1  x  dt 


(A-5) 


Integrating  the  last  term  in  the  integrand  of  Equation  (A-5)  by 
parts,  gives 


I  =  (e[x,t]  -  a'x  } 


t=t. 


t=t  ; tn 

o  o 


{ H [x ,u ,A ,t]  +  A '  x}dt 


(A-6) 


For  small  variations  in  the  control  vectors  (and  consequently 
the  state  vectors)  from  their  optimal  values,  the  resulting  devia¬ 
tion  in  I  is  given  by 


Al  =  I[x  +  6X,  X  +  6x,  U  +  Au]  -  I  [x  ,X  ,U_] 


=  {e[_X  +  6X_ ,t]-A_‘  (x_  +  6X_)  } 


t=t. 


t=t 


{ H [x  +  6X ,  X  +  6X , 


0 


U_  +  6U_, t]  +  X'  (x  +  6x)}dt  -  {e[x,t]  -  A_*  x} 


{H[x,x,u_,t]  +  V  X}dt 


t=t. 


t=t 


0 


(A- 7) 


A-3 


A  Taylor  series  expansion  is  then  made  and  truncated  after  the  first 
order  terms(i.e.  6x_,  6uJ .  The  resulting  change  in  I  is  referred  to 
as  the  first  variation  of  I,  denoted  by  61,  and  is  given  by 


61  =  {6x'[|i  -  xj} 


ax 


t=t 


t=t 


■  aH 


+2j  +  6u'[|]}dt 
o 


(A-8) 


A  necessary  condition  for  an  extremum  of  I  is  that  61 
vanish  for  arbitrary  6_x  and  6u.  Thus  it  is  obtained 


x 


3H 

ax 


(A-9) 


x  =  f(x,u,t)  = 


aH 

au 


0 


and  the  transversal i ty  conditions 

<$x'  [||  -  A]  =  0  at  t  =  tQ,  tf 


(A- 10) 

(A-l 1 ) 


(A-12) 


Estimation  Problem 

The  Bolza  problem  in  the  calculus  of  variations  has  just  been 
discussed.  Now  consider  a  nonlinear  system  with  the  dynamic  equation 

x(t)  =  f(x,t)  +  K(x,t)  w(t)  (A-l 3) 

and  the  observation  equation 

yjt)  =  +  v(t) 


(A-l 4) 


. 
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where 

x_  =  n-state  vector  (including  unknown  parameters) 
f(.x,t)  =  n  vector  function  ,  h_(x,t)  =  m  vector  function 

J<(x., t)  =  n  x  l  matrix  ,  y_(t)  =  m  vector  of  outputs 

w_( t )  =  l  vector  of  unknown  v_(t)  =  m  vector  of  measurement 

inputs  ,  errors 

and  the  control  vector  is  included  in  f_(x_,t). 

It  is  desired  to  find  the  least  squares  estimate  of  x(t), 

A 

designated  by  x_(t),  such  that  the  following  cost  function  is  mini¬ 
mized. 


I 


ftf 

till  -  h(i.t)lln(t)  +  Hi  -  Ki’Oll^t)]  dt 

n  —  ~ 


f  f 


2  2 

[111-  h(i.t)lla(t)  +  Hvtx.t) 


( A- 1 5 ) 


where 

V(x,t)  =  K' (x,t)  W(t)  K(x,t)  (A-16) 

and 

||-  ||^  is  the  "Euclidean  Norm". 

Then  for  the  cost  function  (A-15)  and  system  Equation  (A-13), 
the  Hamiltonian  is 


H(x,A,t)  =||y.(t)  -  h(x>t)||g(t)+||w(t)|||^jt)+X.'(t)[l(x>t)+K(x,t)w(t)] 


(A- 17) 


c 


A-5 


From  the  previous  discussion,  by  setting  ^-=0,  it  follows  that 

oW  — 


2V(x,t)  w(t)  +  K'  (x,t)  =  0  ( A- 18) 

Assume  V_  is  not  singular,  then 

w(t)  =  -}|1(x,t)|,(x,t)i  ( A- 19) 

Substitute  Equation  (A-19)  into  Equation  (A-17) 

H(x,A,t)  =  ||^(t)-h(x,t)||q(t)+A'(t)f(x,t)-i-  A'CtjKCx.tJV'hx.t) 

A_(t)  (A-20) 

Then  as  discussed  before,  the  necessary  conditions  to  have  the  cost 
function  in  Equation  (A-15)  obtain  a  minimum  value 

0  _  3H(j,l,t) 

-  d\_ 

’  =  3H(x,l,t) 

■  A 

3X_ 

Since  tf  is  fixed  and  x(0),  x(tf)  are  free,  so  the  transversal i ty 
condition,  Equation  (A-12),  yields 

A (0 )  =  0  ( A-23) 


(A-21) 

(A-22) 


A_(  t^) 


0 


(A-24) 


, 
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With  Equations  (A-20)  to  (A-24),  a  two-point  boundary  value 
problem  is  formulated.  To  solve  this  kind  of  problem,  Detchmendy 
and  Sridhar  [22]  used  the  invariant  imbedding  technique  and  replaced 
the  boundary  conditions  in  (A-26)  and  (A-27)  by  the  more  general 
conditions 


1(0}  =  0  (A-25) 

X(tf)  =  C  (A-26) 

Let  the  missing  final  conditions  be  given  by 

x(tf)  =  r(C,tf)  (A-27) 

with  1  and  t^  considered  to  be  independent  variables.  Then  by 
Taylor's  expansion,  the  left  hand  side  of  Equation  (A-27)  becomes, 

x_(tf+Atf)  =  x(tf)  +  f(x,t)At  +  A^  (A-28) 

and  the  right  hand  side  of  (A-27)  becomes, 

3  r  3  r 

r(C+AC,tf+Atf)  =  r(C,tf)  +  AC  +  ^  Atf  +  A \  (A-29) 

and  from  Equation  (A-26), 

AC  -  -  f  Atf  +  A23  (A-30) 

where  a^ ,  A^  and  a^  represent  higher  order  terms  of  At^,  At^  and  a£, 
and  a t ^  respectively  and  approach  zero  by  taking  limits. 

Substituting  Equation  (A-30)  into  Equation  (A-29),  equating  the 
right  hand  sides  of  Equation  (A-28)  and  Equation  (A-29),  and  taking 


c 


' 
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limits  gives  the  following  invariant  imbedding  equation 

3r  dr  3H(r,C,tf)  3HCrsC,tf) 

3tf  '  [3?  - W -  =  - (*-31) 

Substituting  Equation  (A-20)  into  Equation  (A-31)  gives 
dr_  dr 

atf  ‘  ^{"2  [y.(tf]  -  h(r,tf)]  +  f'r(r,tf)  C 

-  |c'  ~  [|(r,tf)  V-1(r,tf)  K'  (r,tf)]  C} 

=  f(r,tf)  -  1  K(r,tf)  V_1(r,tf)  K'(r,tf)  C  (A-32) 

where  subscripts  represent  taking  derivatives  of  respective  functions. 

T22l 

Assume  an  approximate  solution  of  _r(£, t^)  of  the  formL  J 

r(C,tf)  =  x(tf)  -  P(tf)  C  (A-33) 

where  P_(t^)  is  an  n  x  n  matrix. 

By  substituting  (A-33)  into  (A-32)  and  expanding  the  result 
with  Taylor's  expansion  about  r.(0_,t^),  it  is  obtained 

dP ( t^)  dx(t^) 

-  C  +  dt^~  +  £(tf){-  2h'~(x,tf)  £(tf)[y.(tf)  -  h(x,tf)] 

+  2[Mh';(i.tf)  £(tf)(*(tf)  -  h(x.tf))]'  p  £  +  f'j  c  -  [Mr^c)]-p  c 

3X  -  “  32L  _ 

-  J-C1  4  [K(x,tf)  V_1(x,tf)  K'(x.tf)]  £  +  (£'  ^7tK(x.,tf)£  hx.tf) 

^  —  '  —  3X  9X 


'  n  ;  * 


. 
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|'(x,tf)]C)]'P  C}  =  f(x,tf)  -  f'~  £_C  -  I|(x,tf) 

K'(x,tf)  C  +  k^rCKCx.tf)  l\x,  tf)  K(x,tf)  C)]1  PC  (A- 34) 

_  ax  -  i  _  t  _ 

Collecting  terms  of  zero  order  £  and  those  of  order  higher 
than  zero  yields 


dx(t.r) 

— Jjf—  =  f(x,tf)  +  P.(tf){2h'j(x,tf)  S.(tf)[y_(tf)-h(x,tf)]}  (A-35) 


dP(t^)  a 

=  {-LCtflf-j-  f'jKtf)  -  2E(tf)|  [H'x(x» tf)a(tf) 


-1 


(y.Ctf)-h(x,tf))]£(tf)  -  ^i(x>tf)V",(x,tf)K(x,tf) 


+  (terms  containing  £'  or  C_) }  £  (A-36) 

From  Equation  (A-36),  the  terms  in  the  bracket  of  the  right  hand  side 
should  be  equal  to  -dP_(t^.)/dt^.  Since  only  £  =  £  is  of  interest,  the 
filter  equations  are  obtained  as 


f(x,t)  +  P(t  ){2  h'j&t)  £(t)  [y_(t)  -  h(x,t)]} 


dP_(  t ) 
dt 


P(t)f'-Cx.t)  +  f5(x.t)P(t)  +  2  p(t){^rC!L,s(x.t) 
=  x.  “X  —  —  9_x  — 

a(t)(£(t)-h(x,t))]>  +  \  K(x,t)  V_1(x,t)  K'(x,t) 


(A-37) 


(A-38) 
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APPENDIX  B 

In  this  appendix  two  computer  programs  are  listed.  One  is 
the  program  written  to  estimate  the  unknown  parameter  E-j  of  the  single 
reaction  example  in  Chapter  III.  The  other  program  estimates  the 
unknown  parameters  E^  and  E^  of  the  complex  reaction  example  in 
Chapter  IV.  Both  of  the  programs  are  written  for  the  Filter  I 
algorithm.  Modifications  of  the  program  can  be  made  to  investigate 
the  effects  of  noise  level  and  data  sampling  interval  on  filter 
estimation.  A  change  of  one  program  statement  in  the  listed  pro¬ 
gram  results  in  the  Filter  II  algorithm.  For  Filter  III  algorithm, 
a  change  in  a  function  evaluates  subroutine  becomes  necessary  as 
is  indicated  in  the  appropriate  place  in  the  program. 
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C  MAINLINE  SMPX 

C  MAINLINE.  SMPX 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•a*#**  *  *•  *•*•#-»■**-*■«-*■*-  **#***•*-*********•*■#•*-«■*#  #***#•■»(•  ****** 


#  * 

*  MAINLINE  SMPX  * 

■*  - - -  # 


*  * 

*  PURPOSE  * 

*  THIS  PROGRAM  IS  WRITTEN  BY  USING  DETCHMENDY  * 

*  AND  SRIDHAR  NONLINEAR  FILTER  TO  IDENTIFY  THE  * 

*  UNKNOWN  ACTIVATION  ENERGY  OF  A  SINGLE  RE-  * 

*  ACTION  EXAMPLE (CH I EN  AND  ARIS)  * 

•*  * 


*  PARAMETER  DESCRIPTIONS  * 


*  Y  (  1  ) . EXTENT  OF  REACTION  * 

*  Y  (  2  J . STATE  VARIABLE  CORRESPONDING  * 

*  TO  TEMPERATURE  * 

*  Y  (  3  . . STATE  VARIABLE  CORRESPONDING  * 

*  TO  UNKNOWN  El  * 

*  PRMT(l) . LOWERBOUND  OF  TIME  * 

*  PRMT ( 3 ) . ...SAMPLING  INTERVAL  * 

*  PRMT  (2 ) . UPPERBOUND  OF  TIME  * 

*  NDIM . DIMENSION  OF  STATE  VARIABLE  * 

*  HR . FORWARD  FREQ.  FACTOR  * 

*  HB . BACKWARD  FREQ.  FACTOR  * 

*  Q  (  I  *  J  ) . COVARIANCE  DIFFERENTIAL  MATRIX* 

*  OF  NONLINEAR  FILTER  * 

*  CT . THE  ESTIMATE  OF  TEMPERATURE  * 

*  EA . THE  ESTIMATE  OF  ACTIVATION  * 

*  ENERGY  OF  El  * 

*  TM . THE  MEASURED  TEMPERATURE  * 

*  * 


*  SUBROUTINE  CALLED  AND  USAGE 

*  CALL  RKGM 

*  I.E.  RUNG E-KUTTA “GILL  INTEGRATION  ROUTINE 


*  CALL  GAUSS ( IX*S*AM*V)  * 

*  I.E.  GAUSSIAN  RANDOM  NUMBER  GENERATOR  * 

*  (IBM  360  SSP  ROUTINE)  5 


*  W  H  E  R  E  ’c 

*  IX  IS  AN  ARBITRARY  I NTEGER ( POSI T I VE )  * 

*  S  IS  THE  STANDARD  DEVIATION  * 

*  AM  IS  THE  AR I  THEMATIC  MEAN  * 

*  V  IS  THE  GENERATED  RANDOM  NUMBER  * 

*  * 

*  ##  *  ****■}<■  #  -«■  **  #  *  **  *■*  ******  ****************  **** 


INTEGER  PA 

DIMENSION  RRX( 1700) *RCA( 1700) »RTP( 1700) 
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C  MAINLINE  SMPX  ...(CONT'D) 

COMMON  PRMT (3) »X*NDIM»Y(12) *DERY ( 12 ) 

COMMON  R  » HR  tHB  » Q ( 3  »  3 )  *Cj*CX»ZH(3)  #CCX*TM#QINT 


C  SPECIFY  THE  VALUES  OF  PRMT ( 1 ) > PRMT ( 2 ) » PRMT ( 3 ) 


PRMT ( 1 ) =0 • 
PRMT ( 2 ) =1700 • 
PRMT ( 3 ) = 1 • 

X  =  PRMT  (  1  ) 


C  SPECIFY  THE  INITIAL  STATE  VARIABLES 


Y (  1  )  =0 • 

Y  (  2  )  =0  • 

Y ( 3 ) =0  •  1 
WRITE (6 #11) 

11  FORMAT (  • 0 • *  15X  * 'NONLINEAR  FILTER  EST I  AMT  ION  OF  STATES 

*  AND  UNKNOWN 
^PARAMETER ' ) 

WR I TE ( 6  *  2  2 ) 

2  2  FORMAT  (///  »3X  *  '^PROCESS'  » 1 7X  »  ' SET* '  *  4X  >  '-^ESTIMATION' 

*  *  2  7X  > 'SET*' /'0 

*  '  *  5X ♦  'TIME'  »  5 X  » 'CONC. '  »7X» ' TEMP • '  » 6X *  ' T I  ME ' * 6X * 'CONC.  ' 

*  >7X*  ' TEMP.  '  » 

*9X#'Q(1#1>  '  *5X*'Q(2*2)  '  >  4X  » ' Q ( 3  »  3 )  '/) 

R  =  2  . 

HR  =  EXP ( 12. ) 

HB=EXP ( 30 • ) 

C J-400 • /  100 • 


C  SPECIFY  THE  INITIAL  COVARIANCE  DIFFERENTIAL  MATRIX 


DO  9  1=1*3 
DO  9  J=  1  *  3 
Q  (  I  *  J )  =  0 • 

9  CONTINUE 

Q( 1*1 ) =  0 • 00  2  5 
Q(2*2)=0.01 
Q  (  3  *  3  )  =  1  • 


C 


CHANGE  THE  Q(I»J)  MATRIX  TO  STATt  VARIABLES  OF  FILTER 


(C’T  .3) ... 

(  X  )  Y  IC  *  *  (  €  J  7  *y 3  ^0MVG} 
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•  S*fl 
{ •SI)qx3*qH 
(  *C£  )  X? *-5H 
.GOI\.COA*lD 


ia  33..AlflAvC3  JAITI/U  3HT  Y3 1 33q3 


£•1*1  ?  oa. 

£fl*v  -  ;0 

•  0  *  (  L  t  I  J  0 
J'J/IT/103  9 
. .  (1*1 )  n 

I0*U«(S#S)0 


T  XMTAM  (Lt  I  )C  3-»T  3DHAH3 


c 


mainline  smpx 


. . . (CONT'D) 


<  =  0 

DO  10  1=1,3 
K  =  K  +  2 

DO  10  J=  1 ,3 
I  JK= I+j  +  K 
Y (  I  JK ) =Q ( I »  j ) 
10  CONTINUE 


SPECIFY  the  parameter  values  of  subroutine  gauss 


I  X= 1 1 1 1 
S=  1  • 

AM  =  0  • 

PA=PRMT ( 2 ) /PRMT { 3 ) 
N=  PA 

X  =  PRMT ( 1 ) 

ZH ( 1 ) =0  • 

ZH ( 2 ) =0  • 

ZH ( 3 ) =0  . 

EX  =  5  • 

CX  =  5  • 


READ  IN  ALL  THE  PROCESS  DATA 


I  J  =  1 

40  K J  = I J  +  4 

READ( 5 ,33 )  (RRX ( I )  >RCA( I )  , RTP ( I ) , I  =  I J , K J ) 
33  FORMA T ( 3 ( 3 A4 )  ) 

I  J  =  KJ+1 

IFUJ.EQ.  1700)60  TO  70 
GO  TO  40 
70  INT=PRMT(3) 

MN  =  Q 


READ  IN  THE  UPDATE  INTERVAL ( IQ  IN ) FOR  ThE  FILTER  2 
ALGORITHM 


IF  I Q I N  =  P R M T ( 3 )  MEANS  FILTER  1  ALGORITHM 


READ( 5 ,71  ) IQIN 
71  FORMAT (14) 


{  : »  t  : ) . . . 
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I X  1 1 *XI 
•  W 

•  0*MA 

l  £  )  T  5  \  (  $  )  TM*i<  *Aq 

1 1 )  t  .'pq«x 

•0«(XIHS 

•  0*<5JHS 
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■  -  .  T  '  A 
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C  MAINLINE  SMPX  .  .  .  (  COIMT  '  D  ) 

ISET=IQIN 
60  MN*MN+ I  NT 
X  A  =  R  R  X  {  M,N  ) 

CA=RCA ( MN ) 

TP=RTP ( MN ) 

CALL  GAUSS (  IX*S»AM»V) 

T  M  =  T  P  +  V 

IF ( MN • EQ • IQIN)NDIM=12 
IF (NDIM.EQ. 12 ) I Q I N= IQ  I N+ 1  SET 
IF (NDIM.EQ12 ) GO  TO  65 
N  D  I M  =  3 
65  CONTINUE 
CALL  RKGM 

IF (X.LT.EX )GO  TO  60 
EX=EX+5 . 

CT=50Q.+100.*Y ( 2 ) 

EA=  (  1  •  +Y  (  3  )  )  #24000  • 

WRITE (6*43) XA*CA> TP *XtY(l)#CT*EA»Y(4).Y(a)*Yll2) 

43  FORMAT (  1  ‘  * 3X * F7 . 2  * 3X * Fb . 4 » 3X » F9 . 3  * 3X * F7 . 2 » 3X * F8 . 4  * 3X 

*  >F9*  3  »  3X  *E13* 

*5»3(2X»E10.3) ) 

I F ( X • GE • PRMT ( 2 ) )GO  TO  100 
IF (CT.LT.300 )GQ  TO  100 
GO  TO  60 
100  STOP 
END 


( v  '  T  D) ... 


< S) Y*.OOI+ *00e«TD 
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SUBROUTINE  FCT 


C 


SUBROUTINE  FCT 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 
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* 

■K 


SUBROUTINE  FCT 


* 


* 


* 


*  PURPOSE 


* 


* 


THIS  SUBROUTINE  INCLUDE  THE  NONLINEAR  FILTER  * 


SUBROUTINE  FCT 

DIMENSION  F  X ( 3 )  *DFX(3,3)  » FQ ( 3 > 3 )  * QF ( 3 ♦ 3 )  » ZC ( 3  *  3  )  » QZ ( 3 
*  *3  )  *  GZQ (3*3) 

DIMENSION  EQ(3>3) *EY(12) *QH(3) »EFX(2) 

COMMON  PRMT ( 3 )  *  X  » ND I M * Y ( 1 2 ) *  DER Y (12) 

COMMON  R  * HR ♦ HB  * Q ( 3  *  3 ) 

COMMON  Cj*CX*ZH(3) *CCX*TM 
CR=(3.-Y(  1)  )*(2.-Y(l)  ) 

CB= ( 1 .+Y ( 1 ) ) **2 
CFR=-5 • +2 • *Y ( 1 ) 

CFB  =  2 • * ( l.+Y( 1 )  ) 

EP=24000. 

CT=50Q.+1Q0.*Y ( 2 ) 

EA= ( 1 • +  Y ( 3  )  )*EP 
ACR=EXP ( ( -EA/CT ) / R ) 

BE=EA+24000 • 

ACB=EXP ( (-BE/CT ) /R) 

FX (  1  ) =HR*ACR#CR“HB*ACB*CB 
FX(2)=CJ*FX(1) 

FX ( 3 ) =0. 

TE=500 • +100 • #Y ( 2 ) 

ZH(2) = ( TM-TE ) / 100 • 

IF ( ND IM*EQ. 12 )G0  TO  110 


C  IF  ND  IM*  1'2  MEANS  FILTER  1  ALGORITHM 


FILTER  3  ALGORITHM  FOR  EVLUATION  OF  Q(I*J)  IS  PUT  IN 
HERE 

DO  120  1=1*3 
120  QH ( I ) =0 • 

DO  125  1=1*3 


I J.  iT  0,-013  ITU  J  t  IHT 

I  -  i  T  1 
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C  SUBROUTINE  FCT  ...(CONT'D) 

DO  125  J= 1 »  3 

125  QH ( I ) =QH ( I )+Q( I  * j )*ZH( J) 

DO  130  1=1*3 
130  DERY ( I )=FX( I )+QH( I ) 

GO  TO  99 
110  RT2=R*CT**2 
EER=EA/RT2 
EEB=BE/RT2 
RT 1  =  -R*C  T 


C  FOR  FILTER  3  ALGORITHM  THE  FOLLOWING  IS  DELETED 


C  CALCULATE  THE  JACOBIAN 


DFX ( 1.1 ) =HR*ACR*CFR-HB*ACB*CFB 

DFX (1*2) = (HR*ACR*CR*EER-HB*ACB*EEB*C8)*100. 

DFX ( 1*3  )  =  (  (HR*ACR*CR) /RT1-(H6*ACB*CB) /RT1 )*( 24000. ) 
DFX ( 2  *  1 )=CJ*DFX(1*1) 

DFX(2*2 )=CJ*DFX( 1 *2) 

DFX ( 2  > 3 ) =C J*DFX (1*3) 

DFX ( 3  *  1 ) =0 • 

DFX ( 3 . 2 ) =0 • 

DFX ( 3  *3 ) =0. 

K  =  0 

DO  100  1=1.3 
K  =  K  +  2 

DO  100  J= 1  *  3 
I  JK= I + J  +  K 
Q(  I  * J  )  =Y (  I  JK) 

100  CONTINUE 

DO  10  1=1.3 
DO  10  J= 1 *3 
FQ (  I  .  J)=0. 

GF ( I »  J  )  =0 • 

10  CONTINUE 
DO  20  K= 1  *  3 
DO  20  1=1.3 
DO  20  J=1 .3 

FQ ( K . I  ) =FQ ( K  » I ) +DFX (K . J  J  *Q ( j . I ) 

20  CONTINUE 
DO  30  K= 1 »  3 
DO  30  1=1.3 
DO  30  J= 1  * 3 

QF ( K . I  )  =  QF ( K . I )+Q(K»j)*DFX( I.J) 

30  CONTINUE 
DO  40  1=1*3 
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SUBROUTINE  FCT 


. . . (CONT ' D ) 


OH ( I ) =0 . 

DO  40  J  =  1 »  3 
ZC ( I  * J  )=0. 

QZ ( I  * J ) =0. 

QZQ ( I » J ) =0. 

40  CONTINUE 
ZC ( 2  *  2 )  =  1  • 

DO  50  K=  1  *  3 
DO  50  1=1*3 
DO  50  J= 1  *  3 

QZ(K* I )=GZ(K* I )+Q(K*j)*ZC( J* I ) 

50  CONTINUE 
DO  60  K= 1 »  3 
DO  60  1=1*3 
DO  60  J  =  1  * 3 

QZQ ( K  *  I  ) =QZQ ( K  » I ) +QZ ( K » J ) *Q ( j *  I ) 
60  CONTINUE 
DO  70  1=1*3 
DO  70  j=1 *3 

EQ ( I  *  J ) =FQ ( I »  J ) +QF ( I » j ) -QZQ ( I  *  j ) 
70  CONTINUE 
K  =  Q 

DO  80  1=1*3 
K  =  K+2 

DO  80  J= 1  *  3 
KI J=K+I+J 
E V ( K I J ) =EQ ( I  * J) 

80  CONTINUE 
DO  95  1=1*3 
DO  95  j= 1  *  3 

QH(  I  ) =QH (  I  ) +Q ( I  *J)*ZH( j) 

95  CONTINUE 
DO  85  1=1*3 
EY (  I  ) =FX ( I ) +QH ( I  ) 

85  CONTINUE 

DO  90  1=1*12 
DERY (  I  )  =E Y (  I  ) 

90  CONTINUE 
99  RETURN 
END 


(  ' '  o> . . . 


<•  m  : )  h 

•  Q~  kti)^ 
•  0  *  (  C  • I )OSQ 

.  I  *  C  S  #  S  )  0  S 
£•!*>  08  OC 

e«x*i  oa  oa 
c#x*w  oe  oa 


(  *U  DS”  (  n>  )0+(  i  •>  J  XJ«  {  1 »so 

3u.  IT/OO  Oc 


e«;»i  oa  ca 
t» lsc  Oa  oc 


(I« J,  (C»>)S0+(1*  )0X  * ( 1 1 > )0£0 

3U/ITKOO  oa 


tti»i  ot  oa 
€ti«c  ov  oa 


<u#n  il;  Xutl)  3* { t  *  I  )  C3 

3w  IT/, 00  OV 


€•1*1  08  oa 

£  •  X  *l  8  .a 

*  •  ■  =  '  4  ) 

€•1*1  89  oa 

*  - 

(  c  >HS*  (  L«  I  )  (I)  HC*  ( I  J  HO 

3UZ IT/00  89 
£«X*l  88  OC 
(  I  )HO  <  I  )  X3*  {  I  )  Y3 

jU,  IT/00  88 

sx*x*i  09  oa 

(  I  > Y3»( I )Y*30 


e- 


9 


C  SUBROUTINE  RKGM 

C  SUBROUTINE  RKGM 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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#**###****•**##**•**■#*****#*#**#*#*.**.**  -X  **#****.*#.*•***#■** 


*  SUBROUTINE  RKGM  * 

*  - - — - —  * 

*  * 

*  FUNCTION  OF  THIS  SUBROUTINE  * 

*  * 

*  RUNGE-KUTTA-GILL  INTEGRATION  ROUTINE  * 

*  * 

*  PARAMETER  DESCRIPTIONS  * 

*  . . THE  INDEPENDENT  VARIABLE  * 

*  NDIM . DIMENSION  OF  DEPENDENT  VARI-  * 

*  ABLES  * 

*  Y(I)... . THE  DEPENDENT  VARIABLES  WITH  * 

*  DIMENSION  NDIM  * 

*  PRMT  (  3  ) . INTEGRATION  INTERVAL  * 

*  DERY(I) . DERIATIVE  OF  STATE  FUNCTIONS  * 

*  * 

*  SUBROUTINE  CALLED  AND  USAGE  * 

#  * 

*  CALL  FCT  * 

*  * 

#  * 


*  *  ******  #  *•*  *  *  #  *■**•*  *******  *  *  ***  ********  *  *************  * 


SUBROUTINE  RKGM 


DIMENSION  RJ ( 27 ) * RK ( 2 7 ) » RQ ( 2 7 ) 

COMMON  PRMT ( 3 ) >X  >NDIM* Y ( 12 ) »DERY  (  12  ) 

COMMON  R  *HR  *HB  »Q ( 3  *3 ) *CJ  >CX*ZH ( 3 ) tCCX  »TM  »QINT 
H=  PRMT ( 3  ) 

5  RX  =  X 

DO  10  I  =  1 1  ND I M 
10  Rj ( I ) =Y  (  I  ) 

CALL  FCT 
DO  20  I  =  1 » N  D I  M 
20  RK ( I ) =  H * D E R Y (  I  ) 

DO  30  1  =  1 » N  D I M 
R j { I )=RJ( I ) +0 • 5  *RK ( I ) 

3  0  RQ ( I ) =RK (  I  ) 

X=RX+H/2 • 

DO  40  I  =  1 » ND I M 
40  YU  )  =  R J  (  I  ) 

CALL  FCT 
DO  30  I = 1 » NDIM 
50  RK ( I ) =H*DERY ( I J 
SQ=1. /SORT ( 2.  ) 
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SUBROUTINE  RK.GM 


... ( CONT  *  D ) 


6  0  Rj (  I ) =Rj  (  I ) +  ( 1 • "SQ ) * ( RK ( I  ) -RQ  ( I )  ) 
CA  =  2 • -SORT ( 2 .  ) 

CB  =  -2 . +3 • /SORT ( 2.  ) 

DO  70  I  =  1  »  N  D  I  M 
70  RQ( I ) =CA*RK( I )+CB*RQ( I ) 

DO  80  I = 1 *NDI M 
80  V (  I  )  =RJ  (  I  ) 

CALL  FCT 
DO  90  1=1 »NDIM 
RK  C I ) =H*DERY ( I ) 

90  RJ( I)=RJ( I )+( 1 • +SQ ) * ( RK ( I )-RQ( I ) ) 
CE=2 . +  SQRT ( 2 • ) 

CF=-2.-3./SQRT(2. ) 

DC  100  I  =  1  *  N  D I M 
100  RQl I ) =CE*RK ( I ) +CF*RQ( I ) 

X=  RX  +  H 

DO  110  1  =  1*  IN  DIM 
110  Y(  I  )=RJ  l  I  ) 

CALL  FCT 
DO  120  I  =  1  *  N  D I M 
120  RK ( I )=H#DERY( I ) 

DO  130  I  =  1  *  ND I M 

130  Y (  I ) =RJ ( I ) +RK ( I  )  /6  .  -RQ ( I ) / 3 • 

6  RETURN 
END 
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C  MAINLINE  CMPX 


C  MAINLINE  CMPX 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


**************  -it-  ************  *  ************************** 


*  * 

*  MAINLINE  CMPX  * 

*  -  * 
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PURPOSE  * 

THIS  PROGRAM  IS  WRITTEN  BY  USING  DETCHMENDY  * 
AND  SRIDHAR  NONLINEAR  FILTER  TO  IDENTIFY  THE  * 
UNKNOWN  ACTIVATION  ENERGIES  OF  A  COMPLEX  * 

REACTION  EXAMPLE ( ROSENBROCK  AND  STOREY)  * 


PARAMETER  DESCRIPTIONS 


Y  (  1  ) . CONC.  OF  A  * 

Y  (2  ) . CONC-  OF  C  * 

Y  (3  ) . CONC.  OF  D  * 

Y  ( 4  )  •  * . STATE  VARIABLE  CORRESPOND  I  NG  * 

TO  UNKNOWN  E4  * 

Y(S>) . STATE  VARIABLE  CORRESPOND  I  NG  * 

TO  UNKNOWN  E3  * 

PRMT(l) . LOWERBOUND  OF  TIME  * 

PRMT  (  2  ) . UPPERBQUND  OF  TIME  * 

PRMT  (  3  ) . SAMPLING  INTERVAL  * 

TEMP . REACTION  TEMPERATURE  * 

Q(I,j). . COVARIANCE  DIFFERENTIAL  MATRIC* 

OF  NONLINEAR  FILTER  * 

WF(I*J) . WEIGHTING  MATRIC  OF  MEASURED  * 

QUANTITIES  * 

IQIN . UPDATE  INTERVAL  OF  COVARIANCE  * 

DIFFERENTIAL  EQUATIONS  * 

(Si D I . . ..DIMENSION  OF  STATE  VARIABLES  * 

ACTVE ( I )••••••• ACT  I  VAT  I  ON  ENERGY  * 

FREQ  (  I  ) . FREQUENCY  FACTOR  * 

COFF  *1* . ACTUAL  E  (  I  ) /ASSUMED  E  <  I  > 


SUBROUTINE  CALLED  AND  USAGE 
CALL  RKGM 

I.E.  RUNGE-KUTTA-GILL  INTEGRATION  ROUTINE 


CALL  GAUSS ( I  X  »  SD ♦ AM  *  RDN ) 

I.E.  GAUSSIAN  RANDOM  NUMBER  GtNERATOR 

WHERE 

IX  IS  AN  ARBITRARY  INTEGER ( POSI T I VE )  * 

SD  IS  THE  STANDARD  DEVIATION  * 

AM  IS  THE  ARITHEMATIC  MEAN  * 

V  IS  THE  GENERATED  RANDOM  NUMBER  * 

* 

CONTROL  SCHEME  USED 

IPI  FEEDBACK  CONTROL  SCHEME 
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mainline  cmpx 


•  •  •  ( CONT  •  D ) 


##■&##•■&**•«•**# -,r -jr  **•«■■«•*********■«•*•}<■*■«■********■}♦■***■)(•***#* 


DIMENSION  PV ( 3 )  * Q Y ( 3 ) 

COMMON  X*Y(27)»DERY(27)  »PRMT(3)  »  FREQ ( 5 )  »EK(5)  *Q(6.6) 

*•  *  W  F  (  3  *  3  )  *  HD ( 3  » 

*6 )  *THD ( 6  *3  )  *ZH ( 3  )  *NDIM* ACTVE ( 5 )  »ACTV2 . ACTV4 * ACTV5  >L 

*  *  T  EMP 

COMMON  R C 0 N C ( 3 ) 

READ ( 5  »  5  )  ( ACT VE (  I  )  *  1  =  1  ,3  ) 

5  FORMAT ( 3F5.0 ) 

RE AD ( 3  *  1 0 )  ( FREQ ( I ) *1  =  1,5) 

10  FORMAT ( 5F3. A) 

WR I TE ( 6  *  1 5  ) 

13  FORMAT ( '0 ' *  1QX  * 'ON-LINE  IDENTIFICATION  OF  KINETIC 

*  PARAMETERS') 

WRITE (6  *20 ) 

20  FORMAT ( '0 ' *15X» 'ACTUAL  ACTIVATION  ENERGY  AS 

*  FOLLOWING' ) 

WRITE (6 *25) ( ACT VE ( I ) *1=1*5) 

25  FORMAT ('0«*10X*'E1='*F5.Q»3X»'E2='*F5.0.3X*'E3='»F5.U 
**3X, ' E4= ' *F5. 

*0  *  3X  * ' E  5= • .F5.0 ) 

COFF1=0 .8 
COFF2=l .2 

ACTVE (4 )=C0FF1*ACTVE(4 ) 

ACTVE ( 5 ) =COFF2*ACTVE( 5 ) 

ACTV4=ACTVE (4) 

ACTV5=ACTVE ( 5 ) 

ACTV4  *  ACTV5  ARE  THE  THE  INITIAL  GUESSED  E4  AND  E5 


WR I TE ( 6 ♦ 30  ) 

30  FORMAT {' 0 '* 15X »' MODEL  ACTIVATION  ENERGY  AS  FOLLOWING') 
WRITE (6 *35)  ( ACTVE ( I) *1  =  1.3)  » ACTV4  * ACTV5 
3  5  FORMAT (  'O'  *10X» 1  El* 1 »F6.1*3X» ' E  2  = '  *F6.1*3X» ' E  3  =  * »  F6. 1 

*  *  3X  * ' E4= »  *F6. 

*  1  *  3X  * ' E  5= ' * F6.1  ) 


38 


40 


WR I TE ( 6  *  38  ) 

FORMAT (  'O'  *25X» 'ESTIMATION'  »16X*'SET'  »35X*' PROCESS '  »7X 
*  * 'SET'  ) 

WR I TE ( 6  * 40  ) 

FORMAT  (  'O'  *5X*'TIME'*4X*'C(A)  '  ,4X*'C(C)  '  *4X*'C(D)  '  *9X 


*• ' <2 '  » 1 3X  * ' K4 

* '  , 14X  * 'K5 '  *8X> ' TIME'  >3X.  'C( A)  '  *3X 


'C(C)  '  »3X. 'C(D)  '  ) 


SPECIFY  THE  VALUES  OF  PRMT ( I ) 

PRMT ( 2 ) =1. 

PRMT ( 1 ) =0  . 

PRMT ( 3 ) =0.005 
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C  MAINLINE  CMPX  •••(CONT'D) 

X  =  PRMT  (  1  ) 

C  SPECIFY  THE  INITIAL  EST I  MATS  OF  STATES 

Y  (  1 )  =  1  • 

Y  (  2  )  =0  • 

Y  (  3  5  =0. 

Y  (  4  )  =0  • 

Y ( 5 ) =0  • 

C  SPECIFY  THE  VALUES  OF  COVARIANCE  MATRIX  Q(I*j) 

DO  60  1=1*3 
DO  60  J  =  1  * 3 
Q  (  I  *  J  )  =0 • 5 

60  CONTINUE 
G(4*4)=250. 

0( 5*5 )=250. 

DO  61  1=1*3 
DO  61  J  =  1  *  5 
HD (  I »  J  )  =0  • 

61  CONTINUE 
HD ( 1  *  1 ) = 1 • 

HD ( 2  *  2 )  =  1  • 

HD ( 3  *  3  )  =  1  • 

DO  62  1=1*3 
DO  62  U  =  1  *  5 
THD(U*I )=HD( I »J) 

62  CONTINUE 

C  CHANGE  THE  COVARIENCE  MATRIC  Q(I*j)  INTO  STATE 

C  VARIABLES  OF  FILTER  EQUATIONS 

K  =  0 

DO  65  1=1*5 
K=  K- I +5 
DO  65  J=  I  .5 
I  JK=I+J  +  K 
Y(IjK)=Q(I*J) 

65  CONTINUE 

C  WEIGHTING  FACTOR 

DO  70  1=1*3 
DO  70  J  =  1  *  3 
WF ( I  *  J ) =0  • 

70  CONTINUE 
WF ( 1  *  1 )  =  1  • 

WF ( 2  *  2 )  =  1  • 

WF ( 3  *  3 )  =  1 • 
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C  MAINLINE  CMPX  . ..(CONT'D) 

DO  75  1=1*3 
ZH ( I ) =0  . 

75  CONTINUE 
TEMP=1120. 

I  X  =  9 
SD  =  0 •  1 
AM  =  0  • 

SPECIFY  THE  UPDATE  INTERVAL ( IQIN )  OF  FILTER  2 
ALGORITHM 

IF  I Q I N  = 1  MEAS  FILTER  1  ALGORITHM 

I Q I N=  1 
I  S  E  T  =  I  Q  I  N 
I T I ME=0 
PX  =  X 

DO  76  1=1*3 
PY ( I ) =Y  (  I  ) 

76  CONTINUE 
74  CONTINUE 

ITIME=ITIME+1 

C  L= 1  I.E.  SECTION  TO  GERENATE  PROCESS  DATA 

L=  1 

ND I M=  3 
CALL  RKGM 
QX  =  X 

DO  77  1=1*3 
QY ( I ) =Y  (  I  ) 

77  CONTINUE 
X=PX 

DO  78  1=1*3 
Y  (  I  )  =PY ( I ) 

78  CONTINUE 

DO  79  1=1*3 

CALL  GAUSS( I X  * SD *  AM  * RDN ) 

RCONC ( I ) =QY ( I )*( l.+RDN) 

IF ( ITIME.LT. IQIN ) NDIM=5 

IF ( ITIME.EQ. IQIN )NDIM=20 

IF ( IT IME.EQ. IQIN ) IQIN= IQIN+ISET 


C  L  =  2  I.E.  FILTERING  SECTION 

L  =  2 

79  CONTINUE 
CALL  RKGM 

C  IPI  FEEDBACK  CONTROL  POLICY 


(  '  •  •  • 


t  •  l  =  1  /•-. 

•  0* ( I ) HS  ; 

•  CSXIaC  '3T 

x.oaa 

•  0  3  M  A 

•j  (kiuuj  :  ;  inT 

•  IT  I  POO  JA 

liT  I  Jii JA  I  ;  JIT  jA3  I*.  ICI  31 


I  =.  •  I  01 
/IK  I 3T  32 1 
-033KITI 

£«I3I  dt 

( n  y« ( i ) yq 
i’J/IT/IOD  dT 
aui  iTirtoo  *v 
j>  3K 1 T I  *  3M I T I 

c.  -IX  H  j  f  A^iJ  J.D  OT  >  1T0:.2  .3.1  I*J 


I  3  J 


£«MIGH 
.  ,  A'O 

£tl«l  YY  oa 
(  I  )  Y*=  (  I  )  YO 

£•1*1  8Y  oa 

l  I  )  Y3  * (  1  )  Y 

it.  I  T/ICD  9 V 
•  i*I  e-T  C  2 

• 

l/3flsl)*II)  YO*  t  1  J  3/OOP 
<  I  )I.TJ.i^ITl ) 3 1 
*  i  .  (  .  .u.  r  i  i3i 

i2  K«1  *;AI01  <  /.  1 0 1  .03. 3N*  I  T  I  )  3 1 


_ 

ITO.re  £  l  TJ11  .3.1  S*J 


if  IT/lOO  S>Y 
JJAD 


Y01J03  J0>  T /  O0  >0A9Q333  131 


c 


MAINLINE  CMPX 


(CONT'D) 


TEMPI = TEMP 

YB=1.-Y (1)-Y(2)-Y(3) 

YC  =  Y ( 2  ) 

A=FREQ (4 ) *ACTVE (4  ) *YB 
B  =  FREQ ( 5 ) *ACTVE ( 5 ) *YC 
AB=A/B 

IF ( AB.LT.O. )GO  TO  31 

I  F (  IDEN  •  EQ#  2 ) GO  TO  22 

IF (X.EQ.O. )GO  TO  12 

PA=(ACTVE (4)-ACTVE(5) ) /ALOG(AB) 

PB=l./PA+l./658. 

TEMP= 1 • /PB 

I F ( TEMP •  GE • 1 1 20 • ) GO  TO  12 
IF ( TEMP.LT.O. ) GO  TO  12 
I  DEN  =  2 
GO  TO  32 
12  CONTINUE 
TEMP=1120. 

GO  TO  32 
22  CONTINUE 

PA=(ACTVE (4)-ACTVE(5)  ) / ALOG ( AB  ) 
PB  =  1 • /PA+ 1 • / 6 58  • 

TEMP= 1 • /PB 

IF ( TEMP.LT.850. ) TEMP=850. 

GO  TO  32 

31  CONTINUE 
TEMP=TEMP1 

32  CONTINUE 

IF ( TEMP. GT. TEMPI ) TEMP= TEMPI 
IF (TEMP.LT.850. ) TEMP=850. 

PX  =  X 

DO  80  1=1*3 
PY( I )=Y( I ) 

80  CONTINUE 
X  =  QX 

DO  81  1=1*3 
Y (  I ) =QY  (  I  ) 

81  CONTINUE 

I F ( X • LT • PRMT ( 2 ) )G0  TO  74 
100  STOP 
END 
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C  SUBROUTINE  FCT 

C  SUBROUTINE  FCT 

Q  •*  ■*  -ft  *  *  *****  K  -X  *•>(•*  -Jf  *  *•  ft  *  ft  ft  -j-r  ft  *  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft-  ft  ft  ft  ft-  ft-  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  *  *  -it-  ft 


C  *  * 

C  *  SUBROUTINE  FCT  * 

C  *  * 

C  *  * 

C  *  PURPOSE  * 

C  *  THIS  SUBROUTINE  INCLUDE  THE  NONLINEAR  FLTER  * 

C  *  EQUATIONS  AND  SUPPLY  THE  REQUIRED  NUMERICAL  * 

C  *  DATA  FOR  SUBROUTINE  RKGM  FOR  INTEGRATION  * 

C  *  ft 


q  *ft*ft*ftft#***ft***ft*ft****ftft*ft*ftft**ft**ftft**ftft*ft**ft*ftftft*ft#** 

SUBROUTINE  FCT 

DIMENSION  WZ ( 3 )  *  GH ( 6  »  3 )  *  QW ( 6 )  *  E  Y (27)  » F  X ( 6 ) »DFX(6>6) 

*  *  FQ ( 6 ♦ 6 ) » QF ( 6  » 6 

* ) *QWH(6*6) *  Q  WQ ( 6  >  6  )  >EQ(6»6)  *  W  H ( 3  » 6 ) 

COMMON  X  » Y ( 27 ) *DERY ( 27 )  >PRMT(3) *  FREQ (5)  *EK(5)  ,Q(6*6) 

*  » W F ( 3  *  3 ) ♦ HD ( 3  » 

*6)  ♦ THD ( 6  *  3  )  >  ZH ( 3  )  »NDIM*ACTVE(5)  > ACT V2 *  AC T V4 *  ACT V5 » L 
* ♦ T  EMP 

COMMON  RCONC ( 3 ) 

CONST =1 . /TEMP-1 . /658. 

EK(1)=FREQ(1)*EXP( -8000.#CONST ) 

EK ( 2 ) =FREQ ( 2 ) *EXP ( -700Q.*CONST ) 

EK ( 3  > =  FREQ (  3 ) *EXP ( -7500. *CONST ) 

GO  TO  (1*3)  *  L 

1  CONTINUE 

EK (4 ) =FREQ (4 ) *EXP ( -5000. *CONST ) 

EK (5 ) =FREQ ( 5 ) ftEXP ( -7500 . *CQNST ) 

FX ( 1 ) =- ( EK ( 1 ) +EK ( 2 ) +  EK ( 3 ) ) *Y ( 1 ) 
FX(2)=(1.-(Y(1)+Y(2)+Y(3) ) )*EK(4)-EK(5)*Y(2) 

FX ( 3 ) = ( EK (2 ) +EK ( 3 ) ) *Y ( 1 ) +  EK ( 5 ) *Y ( 2 ) 

DO  2  1=1*3 
DERY ( I  )  =FX  (  I  ) 

2  CONTINUE 
GO  TO  141 

3  CONTINUE 

ACTVE  (  4 ) =ACTV4*Y ( 4 ) +ACTV4 
ACTVE ( 5 ) =ACTV5*Y ( 5 ) +ACTV5 
EK ( 4 ) =  F  R  E  Q ( 4 ) *EXP ( -ACTVE ( 4 ) *CONST ) 

EK(5)=FREQ(5)*EXP( -ACTVE (5 ) *CONST ) 
FX(1)=-(EK(1)+EK(2)+EK(3) ) *Y ( 1  ) 

FX(2)=(1.-(V( 1)+Y(2)+Y(3) ) ) *EK ( 4 ) -EK ( 5 ) * Y (2) 
FX(3)=(EK(2)+EK(3) )*Y(1) +  EK ( 5 )*Y(2> 

FX ( 4 ) =0  • 

FX ( 5 ) =0  • 

ZH(l)=RCONC(l)-Y(l) 

ZH(2)=RCONC(2)-Y(2) 
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C  SUBROUTINE  FCT  ...(CONT'D) 
ZH(3)=RCONC(3)-Y (3) 

CALCULATION  OF  Q(I.J)  ARE  PUT  IN  HERE  FOR  FILTER  3 
K  =  0 

DO  20  1  =  1.5 
K=K-I+5 
I  J=  I 

DO  20  j  = I J  *  3 
I  JK= I +J  +  K 
Q (  I  *  J )  =  Y (  I JK) 

20  CONTINUE 
DO  30  1=1*5 
I  J=  I 

DO  30  J  =  I  J  *  3 
Q( J* I ) =Q ( I > J ) 

30  CONTINUE 
DO  40  J= 1 *3 
WZ ( J) =0  . 

DO  40  1=1*5 
QW (  I  ) =0  . 

QH ( I  * J ) =0* 

40  CONTINUE 
DO  45  K=  1  ♦  5 
DO  45  1=1*3 
DO  45  j  =  1  »  5 

QH(K*  I  )=QH(K* I )+Q(K*j)*THD(j» I ) 

45  CONTINUE 
DO  50  1=1*3 
DO  50  J=  1 .3 

WZ ( I ) =WZ ( I ) +WF ( I ♦ J )#ZH ( J ) 

50  CONTINUE 
DO  55  1=1.5 
DO  55  J=1 .3 

QW (  I  ) =QW ( I ) +QH (  I  >  J  )  *WZ ( J ) 

55  CONTINUE 

ESTIMATION  EQUATION  OF  STATE  VARIABLE  INCLUDING 

UNKNOWN 

PARAMETERS 

DO  60  1=1.5 

EY ( I ) =  FX (  I  ) +QW (  I  ) 

60  CONTINUE 
IF(NDIM.EQ.20)GO  TO  62 
DO  61  I  =  1 .  ND I M 

DERY (  I  )  =  E Y  (  I  ) 

61  CONTINUE 

IF (NDIM.EQ.5 ) GO  TO  141 

62  CONTINUE 
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C  SUBROUTINE  FCT  ...(CONT'D) 


WHEN  FILTER  3  ALGORITHM  IS  USED  THE  FOLLOWING  IS 
DELETED 

DO  10  1=1*5 
DO  10  J=  1  *  5 
DFX ( I  * j  )  =0. 

10  CONTINUE 

DFX ( 1 *1 )=-(E<(l)+E<(2)+E<(3) ) 

DFX ( 2  *  1 )=-EK(4) 

DFX(2»2)=-E<(4)-E<(5) 

DFX ( 2  *  3 ) =-EK ( 4 ) 

DFX (2*4)= ( 1 • - ( Y ( 1 )+Y(2 ) +Y ( 3 ) ) ) *E< ( 4 ) * ( -ACTV4*C0N$T ) 
DFX (2 *5 )=-Y (2 ) *E< ( 5 ) * ( -ACTV5*CONST ) 

DFX (3*1 ) =E<( 2 ) +EK ( 3 ) 

DFX(3*2)=EK(5) 

DFX ( 3  »  5 ) = Y ( 2 ) *E< ( 5  )  * ( -ACTV5*CONST ) 

DO  65  1=1*5 
DO  65  J=  1  * 5 
FQ( I , j ) =0. 

OF (  I  * J  )  =0. 

QWH ( I  * J ) =0. 

QWQ ( I  * J  )  =0. 

65  CONTINUE 
DO  70  1=1*3 
DO  70  J  =  1 »  5 
WH  (  I  *  J ) =  0  • 

70  CONTINUE 
DO  75  <=1*5 
DO  75  1=1*3 
DO  75  J=  1  *3 

WH ( I  *  < ) =WH (  I  *  < ) +WF ( I  * J ) *HD ( J  *K ) 

75  CONTINUE 
DO  80  <=1*5 
DO  80  1=1*5 
DO  80  J= 1 »  3 

QWH ( K • I ) =QWH ( K  *  I ) +QH ( K  *  J ) *WH ( J  » I ) 

80  CONTINUE 

DO  110  <=1*5 
DO  110  1=1*5 
DO  110  J=  1  *  5 

QWQ l < ♦ I  ) =QWQ ( <  *  I ) +QWH ( <  *  J ) *Q ( J  *  I ) 

110  CONTINUE 

DO  90  <=1*5 
DO  90  1=1*5 
DO  90  J  =  1  *  5 

FQ(<*  I  )=FQ(<» I )+DFX(<* j) *Q( J* I ) 

QF ( <  *  I  )=QF(<*I)+0<<*U)*DFX(I*J) 

90  CONTINUE 
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SUBROUTINE  FCT 


.  .  .  (CONT ' D ) 


ESTIMATION  OF  COVARIANT  MATRIX  ELEMENT 

DO  120  1=1*5 
DO  120  J= 1  *  5 

EQ(  I *J)=FQ(  I *j)+QF( I  * J ) -QWQ ( I  * j ) 

120  CONTINUE 
K  =  0 

DO  130  1=1*5 
K>  =  K-  I  +5 
DO  130  J=  I  »  5 
I JK=K+I +j 
EY ( I jK ) =EQ ( I *J) 

130  CONTINUE 

DO  140  1  =  1*  ND I M 
DERY (  I  )=EY(  I  ) 

140  CONTINUE 

141  CONTINUE 
RETURN 
END 
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C  SUBROUTINE  RKGM 

C  SUBROUTINE  RKGM 


C 

C 

C 
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*********  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ********  *.  *  *  *  *  *********  *  *  *  *  *  * 

*  * 
*  SUBROUTINE  RKGM  * 


*  -  # 

*  * 

*  FUNCTION  OF  THIS  SUBROUTINE  * 

*  * 

*  RUNGE-KUTT  A-G I  LL  INTEGRATION  ROUTINE  * 

*  * 

*  PARAMETER  DECRIPTIQNS  * 

*  NDIM . DIMENSION  OF  DEPENDENT  VARI-  * 

*  ABLES  * 

*  Y(I) . THE  DEPENDENT  VARIABLES  WITH  * 

*  DIMENSION  NDIM  * 

*  X . THE  INDEPENDENT  VARIABLE 

*  PRMT  (  3  ) . INTEGRATION  INTERVAL  * 

*  DERY(I) . DERIVATIVE  OF  STATE  FUNCTIONS  * 

*  * 

*  SUBROUTINE  CALLED  AND  USAGE  * 

*  CALL  FCT  * 

*  CALL  OUTP  * 

*  * 


********  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ********  *  ********************** 


SUBROUTINE  RKGM 


DIMENSION  R  J ( 2  7  )  »RK(27) *  R  Q  (  2  7  ) 

COMMON  X»Y ( 27 ) »DERY(27 ) * PRMT (3)  »  FREQ ( 5 )  » EK ( 5 ) * Q ( 6  *  6 ) 

*  > WF ( 3  *  3 ) *  HD ( 3  * 

*6 )  >  THD ( 6  *3 ) *ZH ( 3  )  *NDIM  » ACTVE ( 5 ) *ACTV2  t ACTV4# ACTV5  »L 

*  *  TEMP 

COMMON  RCONC ( 3 ) 

H=  PRMT ( 3  ) 

5  RX  =  X 

DO  10  I  =  1 >  N  D I M 
10  RJ ( I ) =Y ( I ) 

CALL  FCT 
DO  2  0  I  =  1 ♦ N D I M 
20  RK ( I ) =H*DERY ( I ) 

DO  30  I  =  1 » ND I M 
R j (  I  )=RJ(  I  )+0.5*RK( I ) 

30  RQ ( I ) =RK ( I ) 

X  =  RX  +  H/2 • 

DO  40  I  =  1  *  ND I M 
40  Y (  I  )  =R J (  I  ) 

CALL  FCT 
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SUBROUTINE  RKGM 


.  .  .  ( CONT • D ) 


DO  50  I  =  1  *  ND I M 
50  RK ( I ) =  H  *  D  E  R  Y (  I ) 

SG  = 1 • /SORT ( 2 • ) 

DO  60  I=1»NDIM 

60  Rj ( I ) =RJ ( I ) + ( l.-SQ)*{RK( I ) -RQ ( I ) ) 
CA  =  2  • -SORT ( 2 • ) 

CB=-2.+3./SQRT(2. ) 

DO  70  I  =  1  *  ND I M 
70  RQ{ I ) =CA*RK( I )+CB*RQ( I ) 

DO  80  I  =  1  *  N  D  I  M 
80  Y (  I  ) =RJ (  I  ) 

CALL  FCT 
DO  90  I=1»NDIM 
RK ( I ) =H*DER Y ( I ) 

90  Rj( I ) =R J ( I )  +  ( l.+SQ)*(RK( I ) -RQ ( I ) ) 
CE  =  2 • +SQR  T ( 2 • ) 

CF=-2.-3./SQRT(2. ) 

DO  100  I  =  1 » N  D I M 
100  RQ ( I )=CE*RK( I )+CF*RQ( I ) 

X=  RX  +  H 

DO  110  I  =  1  *  ND I M 
110  Y (  I  )  =R J ( I ) 

CALL  FCT 
DO  120  I  =  1 f  ND I M 
120  RK (  I  ) =H*DERY ( I ) 

DO  130  I  =  1  *  ND I M 

130  Y ( I ) =R J ( I ) +RK ( I ) / 6 • -RQ ( I  )  / 3 • 

CALL  OUTP 
6  RETURN 
END 
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SUBROUTINE  OUTP 
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SUBROUTINE  OUTP 
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SUBROUTINE  OUTP 


ft 


* 


*  PURPOSE 


ft 


THIS  SUBROUTINE  IS  THE  OUTPUT  ROUTINE  TO 
PRINTE  OUT  THE  REQUIRED  DATA 


* 


SUBROUTINE  OUTP 


COMMON  X»Y(27) »DERY(27) *PRMT(3)  *  FREQ (5)  *EK(5),Q(6»6) 

*  *WF ( 3  *3 )  »HD ( 3  * 

*6)  *  THD ( 6  »  3 )  »  ZH l 3  )  » ND I M * ACTVE ( 3 )  » ACTV2 *  AC T V4 » ACT V5 ♦ L 
* » T  EMP 

COMMON  RCONC ( 3 ) 

GO  TO  ( 1 *2 ) »L 

1  CONTINUE 

WRITE(6*5)X»Y(1)  *  Y ( 2  >  *  Y ( 3 ) *  TEMP 
5  FORMAT ( '  ' *85X»F5.3»3(2X>F5.4) *3X»F5.0) 

GO  TO  11 

2  CONTINUE 
9  CONTINUE 

ACTVE (4 ) =ACTV4*Y (4 )+ACTV4 
ACTVE ( 5 ) =ACTV5*Y ( 5 )+ACTV5 
WRITE ( 6  » 10 ) X  *Y  (  1  )  » Y ( 2  )  *  Y  C  3 ) * ACTVE ( 4 )  » ACTVE ( 5 ) 

10  FORMAT ( '  +  1 ♦3X»F3.3*3X»F3.4»3X*F3.4»3X»F3.4»16X*2(3X 

*  i  E 13 • 5  )  ) 

11  CONTINUE 
RETURN 
END 
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